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1.  SUMMARY 


Recent  advances  in  seismic  theory  have  provided  more  accurate  representation  of  the 
propagation  of  finite-frequency  seismic  waves  than  ray  theory.  In  this  project,  we  use  the 
“banana-doughnut”  sensitivity  kernels  of  teleseismic  body  waves  to  image  the  crust  and 
mantle  beneath  eastern  Eurasia.  We  have  collected  and  processed  available  broadband 
data  from  both  permanent  stations  and  temporary  networks  in  eastern  Eurasia.  In 
southeast  Tibet,  where  a  PASSCAL  experiment  provided  a  dense  station  coverage,  a 
detailed  study  is  carried  out  to  obtain  high-resolution  P-  and  S-velocity  models,  which 
provide  not  only  important  constraints  on  the  continental  collision  and  plateau  building 
processes  but  also  a  reference  for  comparison  with  a  larger  regional-scale  inversion.  The 
P  and  S-wave  velocity  models  reveal  a  low-velocity  anomaly  in  the  crust  and  upper 
mantle  to  ~  300  km  depth  beneath  a  north-south-trending  rift  zone  in  southeastern  Tibet. 
This  low  velocity  anomaly  is  situated  above  a  tabular,  high-dipping-angle,  high-velocity 
anomaly  that  extends  into  the  upper  mantle  transition  zone.  These  results  are  evidence 
for  the  delamination  of  the  mantle  lithosphere  and  its  causal  relationship  to  the  formation 
of  the  north-south  trending  rift  in  southeastern  Tibet.  A  regional  P-wave  velocity  model 
for  eastern  Eurasia  is  also  obtained  by  inverting  the  finite-frequency  traveltimes  of 
teleseismic  body  waves  from  available  broadband  stations  in  eastern  Eurasia.  In 
southeast  Tibet,  the  regional  model  is  quantitatively  consistent  with  the  results  of  the 
local  and  detailed  study.  Because  of  steep  incidence  of  teleseismic  body  waves  at 
shallow  depth,  gaps  in  resolution  exist  in  the  shallow  upper  mantle  and  crust  in  places 
with  sparse  stations. 

In  the  crust  and  uppermost  mantle,  where  the  velocity  structure  is  highly 
heterogeneous,  a  ID  reference  model  and  single-scattering  approximation  may  become 
invalid.  To  improve  resolution  in  the  crust  and  uppermost  mantle,  we  have  made  several 
contributions  to  the  finite-difference,  scattering-integral  method  (FDSIM).  These  include 
the  finite-frequency  sensitivity  kernels  for  head  waves  and  component-dependent 
sensitivities  due  to  uneven  distribution  of  scattered  waves  on  the  different  components  of 
the  same  arrival.  Furthermore,  we  apply  the  finite-frequency  methodology  to  surface- 
wave  tomography  using  the  ambient  seismic  noise.  Unlike  previous  studies  that  use 
estimated  Green’s  functions  to  solve  for  phase  or  group  velocity  maps  and  then  a  point- 
wise  S  velocity  structure  under  a  ID  assumption,  we  construct  3D  sensitivity  kernels  for 
the  finite-frequency  Green’s  functions,  taking  into  account  a  3D  reference  model  and  the 
effects  of  topographic  variations.  The  sensitivity  kernels  and  the  phase  delays  between 
the  estimated  Green’s  functions  and  synthetics  are  inverted  for  Vp  and  Vs  structure 
beneath  southeast  Tibet.  Results  show  that  the  wave  speed  in  the  reference  model  (a 
modified  CRUST  2.0)  is  too  high  in  the  lower  crust  and  too  low  in  the  shallow  crust 
beneath  southeast  Tibet. 
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2.  INTRODUCTION 


In  seismic  tomographic  studies  to  date  the  reference  models  are  commonly  ID  and 
the  structural  sensitivity  kernels  of  seismic  data  are  calculated  without  considering  the 
finiteness  of  seismic  waves  in  both  time  and  frequency  domains.  These  simplifications 
result  in  a  theoretical  limit  (Nolet  and  Dahlen,  2000;  Baig  et  al.,  2003),  in  addition  to  that 
from  the  data  coverage,  on  the  structural  resolution  that  can  be  achieved  in  tomography. 
In  recent  years,  seismologists  have  started  to  replace  the  great-circle  path  of  a  surface 
wave  with  sensitivity  kernels  having  a  finite  width  and  symmetrically  distributed  across 
the  great-circle  path  (e.g.,  Yoshizawa  and  Kennett,  2002;  Ritzwoller  et  al.,  2002)  or  a 
more  physically  realistic  kernel  calculated  using  ID  (e.g.,  Li  and  Romanowicz,  1996; 
Zhou  et  al.,  2004,  2006)  or  3D  (e.g.,  Liu  and  Tromp,  2006;  Shen  et  al.,  2008b)  reference 
models.  Similarly,  several  studies  (e.g.,  Hung  et  al.,  2004;  Montelli  et  al.,  2004;  2006; 
Yang  et  al.,  2006;  Chen  et  al.,  2007b;  Sigloch  et  al.,  2008;  Ren  and  Shen,  2008)  have 
replaced  body-wave  ray  paths  with  “banana-doughnut”  sensitivity  kernels  calculated  in 
ID  (Dahlen  et  al.,  2000;  Hung  et  al.,  2000;  Zhao  et  al.,  2000;  Zhao  and  Jordan,  2006)  or 
3D  (Tromp  et  al.,  2005;  Zhao  et  al.,  2005;  Liu  and  Tromp,  2006;  Zhang  et  al.,  2007) 
reference  models. 

The  use  of  body  wave  ray  theory,  a  ID  reference  model,  and  Bom  single-scattering 
approximation  makes  it  possible  to  calculate  the  sensitivity  kernels  of  high  frequency 
teleseismic  body  waves  at  a  relatively  low  computational  cost  (Dahlen  et  al.,  2000;  Hung 
et  al.,  2000).  On  the  other  hand,  computationally  intensive  numerical  approaches  such  as 
the  finite-difference  or  spectral-element  method  represent  more  accurately  complex 
interactions  of  full  waves  in  3D  heterogeneous  media  and  make  it  possible  to  iteratively 
improve  upon  a  3D  model. 

In  the  following  sections,  we  document  our  application  of  the  “banana-doughnut” 
sensitivity  kernels  based  on  ID  reference  models  to  teleseismic  body- wave  traveltime 
tomography  in  east  Eurasia.  We  demonstrate  that  high-resolution  tomographic  models  of 
the  crust  and  upper  mantle  can  be  achieved  using  finite-frequency  teleseismic  body 
waves  in  places  with  dense  stations,  such  as  southeast  Tibet.  We  show  that  the 
continental-scale  P  velocity  model  obtained  from  finite-frequency  teleseismic  body  wave 
tomography  is  comparable  to  that  of  the  high-resolution  study  in  southeast  Tibet.  In 
order  to  improve  resolution  in  the  crust  and  uppermost  mantle,  we  made  several 
contributions  to  the  full-wave  tomography  method  based  on  finite-difference  simulations 
of  wave  propagation  in  3D  reference  models.  We  applied  the  full-wave  approach  to 
image  the  crust  and  upper  mantle  beneath  southeast  Tibet  using  ambient  noise. 

3.  TECHNICAL  APPROACH 

3.1  Finite-Frequency  Tomography  of  Teleseismic  Body  Waves 

3.1.1  Finite-Frequency  Traveltime  Data 

We  collected  and  utilized  available  broadband  waveforms  to  construct  the  finite- 
frequency  velocity  model  beneath  eastern  Eurasia.  Figure  1  shows  the  distribution  of  the 


2 


stations,  which  include  the  Global  Seismographic  Network  (GSN),  Japanese  F-net  and 
JISNET,  Taiwan  Broadband  Seismic  Network,  and  other  regional  seismic  networks 
available  at  the  Data  Management  Center  (DMC)  of  the  Incorporated  Research 
Institutions  for  Seismology  (IRIS).  We  selected  global  earthquakes  in  the  updated 
Engdahl-van  der  Hilst-Buland  (EHB)  bulletins  prior  to  2004  and  the  National  Earthquake 
Information  Center  (NEIC)  bulletins  for  recent  events  with  magnitude  greater  than  5.5 
recorded  by  the  seismic  stations  in  eastern  Eurasia.  We  used  event-station  paths  with 
epicentral  distances  between  30°  and  85°.  Delay  times  were  measured  by  an  automated 
waveform  cross-correlation  routine  based  on  VanDecar  and  Crosson  (1990)  in  three 
frequency  bands  for  P  waves  (0. 5-2.0  Hz,  0. 1-0.5  Hz  and  0.03-0. 1Hz)  and  S  waves  (0.1- 
0.2  Hz,  0.05-0.1  Hz  and  0.02-0.05  Hz).  The  signal-to-noise-ratio  threshold  was  20  for 
quality  control  in  each  frequency  band  and  the  selected  records  were  visually  inspected 
for  consistency.  The  signal-to-noise  ratio  is  defined  as  the  ratio  of  the  peak-to-peak 
amplitude  of  the  main  arrival  to  the  standard  deviation  of  the  time  series  in  an  80-s 
window  before  the  main  arrival. 

To  reduce  the  tradeoff  between  crustal  and  mantle  velocity  heterogeneities  in  seismic 
tomography,  we  removed  the  frequency-dependent  crustal  correction  from  the  teleseismic 
travel  times  by  cross-correlating  the  impulse  responses  of  a  crust  model  filtered  in  a 
narrow  frequency  band  (Yang  and  Shen,  2006;  Appendix  I).  We  used  a  crustal  model 
from  Sun  et  al.  (2005)  for  seismic  stations  in  China,  and  CRUST2.0  (Bassin  et  al.,  2000) 
for  stations  elsewhere  to  calculate  the  traveltime  difference  with  respect  to  IASP91. 

3.1.2.  P-  and  S-Velocity  Models  in  Southeast  Tibet 

The  distribution  of  seismic  stations  in  eastern  Eurasia  is  uneven  (Figure  1).  In  places 
where  stations  are  dense,  it  is  possible  to  carry  out  detailed  data  analyses  to  obtain  high- 
resolution  P-  and  S-velocity  models.  In  addition  to  collecting  data  needed  for  the 
continental-scale  inversion,  these  detailed  studies  allow  us  to  establish  references  to 
which  large-scale  models  can  be  compared.  With  several  PASSCAL  experiments  (Figure 
1),  southeastern  Tibet  is  such  a  region  for  a  detailed  study.  High-resolution  images  of  the 
crustal  and  mantle  structure  in  the  region  also  provide  important  constraints  on  the 
continental  collision  and  plateau  building  processes. 

We  processed  data  from  the  Namche  Barwa  Seismic  Experiment  (Sol  et  al.,  2007), 
which  deployed  a  50-station  broadband  array  and  a  20-element  short-period  array  in 
southeastern  Tibet  during  2003-2004.  Data  from  the  GSN  station  LSA  were  also 
included  in  this  study.  Following  the  data  processing  method  outlined  above,  we  obtained 
44,000  P  and  19,500  S  delay  times,  which  were  then  utilized  to  invert  for  spatial 
variations  in  P-  and  S-wave  velocity  perturbations  according  to  the  3-D  finite  frequency 
kernel  formulation  (Dahlen  et  al.,  2000;  Hung  et  al.,  2004;  Yang  et  al.,  2006). 

The  model  space  is  parameterized  with  a  regular  3D  grid  of  33x33x33  centered  at 
(93°  E,  30°  N),  and  has  a  dimension  of  18°  in  longitude,  14°  in  latitude,  and  1200  km  in 
depth.  The  grid  spacing  is  approximately  53  km  in  longitude,  46  km  in  latitude,  and  40 
km  vertically.  The  inverse  problem  is  resolved  using  a  standard  damped-least-square 
algorithm  (Paige  and  Saunders,  1982)  and  the  damping  parameter  is  determined  through 
a  trade-off  analysis  of  model  norm  versus  variance  reduction.  We  use  checkerboard 
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resolution  tests  with  the  same  damping  parameter  as  in  the  real  data  inversion  to  evaluate 
the  ability  of  a  given  data  coverage  and  inversion  technique  to  recover  crustal  and  mantle 
structures.  These  tests  show  that  structures  are  well  recovered  beneath  the  stations  at 
depth  between  -50  km  and  -450  km.  The  minimal  size  of  the  laterally  and  vertically 
resolved  structure  is  ~  100  km. 

Our  results  (Ren  and  Shen,  2008;  Appendix  II)  show  a  prominent  N-S  trending  low- 
velocity  anomaly  at  ~92°E  down  to  at  least  300  km  depth  (Figure  2).  These  structures 
are  observed  on  both  Vp  and  Vs  models.  We  associate  the  high-velocity  anomalies 
southwest  and  east  of  the  north-south  trending  low  velocity  anomaly  to  the  Indian 
lithosphere,  which  extends  to  -250  km  depth.  We  observe  no  evidence  for  a  subducted 
Indian  lithosphere  at  depth  of  300  km  or  more  beneath  the  study  area.  The  north-south 
low  velocity  anomaly  across  Indus- Yarlung  Suture  coincides  strikingly  well  with  a  rift  on 
the  surface  at  southeastern  Tibet  (e.g.,  Tappnnier  and  Monlar,  1977;  Armijo  et  al.,  1986; 
Yin,  2000)  and  is  above  a  high-angle,  southeastward  dipping,  high-velocity  anomaly  that 
extends  into  the  upper  mantle  transition  zone,  as  shown  on  the  cross-sections  (Figure  2c). 
We  associate  the  high-angle  dipping,  high-velocity  anomaly  to  a  sunken,  delaminated 
mantle  lithosphere  and  the  north-south-trending  low-velocity  structure  asthenospheric 
upwelling  induced  by  the  sunken  mantle  lithosphere.  The  upwelling  asthenospheric 
mantle  may  be  responsible  for  the  thermal  weakening  of  the  overlying  crust  and  the 
localization  the  rifts  observed  on  the  surface  of  the  Tibetan  Plateau.  Our  velocity  models 
suggest  that  the  north-south  trending  rifts  on  the  Tibetan  plateau  are  lithosphere-scale 
features  and  the  east-west  extension  of  the  plateau  is  mechanically  coupled  throughout 
the  crust  and  mantle  lithosphere.  The  relation  between  north-south-trending  rifts  on  the 
plateau  and  underlying  low-velocity  anomalies  is  confirmed  elsewhere  in  southern  Tibet 
(Liang  et  al.,  2008). 

3.1.3.  P-velocity  model  of  East  Eurasia 

We  also  carried  out  inversions  for  the  P  velocity  structure  beneath  eastern  Eurasia. 
The  tomographic  model  extends  to  2500  km  depth  and  is  parameterized  by  a  81x81x51 
grid  with  spacing  of  1.0°,  1.25°,  and  58  km  in  latitude,  longitude  and  depth,  respectively. 
As  in  the  detailed  study  in  southeast  Tibet  (Ren  and  Shen,  2008),  this  inversion  is  carried 
out  using  the  finite-frequency  tomography  methodology  (Hung  et  al.,  2004;  Yang  et  al. 
2006).  We  also  utilize  a  convolutional  quelling  technique  (Meyerholtz  et  al.,  1989)  in  this 
continental-scale  study. 

Figure  3  shows  the  P-velocity  anomalies  at  180,  270,  and  540  km  depth  beneath  east 
Eurasia.  The  model  in  southeast  Tibet  is  generally  consistent  and  comparable  to  the 
images  from  the  detailed  study  (e.g.,  the  north- south- trending  low  velocity  anomaly). 

An  inspection  of  the  sampling  of  the  structure  beneath  eastern  Eurasia  by  the  finite- 
frequency  sensitivity  kernels  reveals  that  there  is  a  good  data  coverage  in  the  depth  range 
of  350  -  1200  km  in  most  parts  of  the  study  area,  but  gaps  exist  in  the  shallow  upper 
mantle  and  crust,  except  in  regions  with  dense  stations  (e.g.,  southeast  Tibet,  Japan). 
Models  based  on  teleseismic  P  waves  have  limited  vertical  resolution  in  the  crust  and  uppermost 
mantle,  because  of  steep  incidence  angles  and  lack  of  crossing  rays  at  shallow  depths.  The 
problem  is  particularly  severe  in  places  of  sparse  station.  To  obtain  a  good  resolution  at  shallow 
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depth,  we  need  to  incorporate  surface  waves  and  regional  body  waves.  In  the  crust  and 
uppermost  mantle,  however,  the  velocity  structure  is  highly  heterogeneous  and  the  ID 
reference  model  used  in  the  calculation  of  the  “banana-doughnut”  sensitivity  kernels  may 
become  invalid.  It  is  thus  necessary  to  go  beyond  the  ID  reference  model  and  single¬ 
scattering  approximation  (Dahlen  et  al.,  2000)  and  to  develop  a  full- wave  approach  that 
accounts  for  complex  wave  interaction  in  the  crust  and  shallow  mantle.  The  following 
sections  report  the  development  of  the  finite-difference  scattering-integral  method 
(FDSIM)  and  its  application  to  image  the  earth  structure  beneath  eastern  Eurasia. 

3.2.  Finite-Frequency  Sensitivity  Kernels  for  Pn/Sn  Waves 

We  examined  the  finite-frequency  effects  and  sensitivity  kernels  of  head  waves 
(Zhang  et  al.,  2007),  which  are  extremely  important  in  determining  the  structure  of  the 
predominantly  layered  Earth.  A  numerical  experiment  was  designed  to  demonstrate  the 
finite-frequency  effects.  The  model  has  a  low-velocity  layer  over  a  high-velocity  half 
space  and  a  cylindrical-shaped  velocity  anomaly  placed  beneath  the  interface  at  different 
locations.  A  3D  finite-difference  method  (Olsen,  1994)  is  used  to  calculate  synthetic 
waveforms.  Travel-time  and  amplitude  anomalies  are  measured  by  the  cross-correlation 
of  synthetic  seismograms  with  and  without  the  velocity  perturbation  and  are  compared  to 
the  3D  sensitivity  kernels  constructed  from  full  waveform  simulations.  The  results  show 
that  the  head  wave  arrival-time  and  amplitude  are  influenced  by  the  velocity  structure 
surrounding  the  ray  path  in  a  pattern  that  is  consistent  with  the  Fresnel  zones.  Unlike  the 
“banana-doughnut”  travel-time  sensitivity  kernels  of  turning  waves,  the  travel-time 
sensitivity  of  the  head  wave  along  the  ray  path  below  the  interface  is  weak,  but  non-zero 
(Figure  4).  Below  the  ray  path,  the  travel-time  sensitivity  decreases  with  depth  from  the 
interface  to  a  minimum  (zero)  before  increasing  to  a  maximum  (absolute  value),  the 
depth  of  which  depends  on  the  wavelength  and  propagation  distance.  Thus  head  waves 
with  horizontally  crossing  ray  paths  do  not  necessarily  sample  the  same  structure  beneath 
the  interface.  The  sensitivity  kernels  vary  with  the  vertical  velocity  gradient  in  the  lower 
layer,  but  the  variation  is  relatively  small  at  short  propagation  distances  when  the  vertical 
velocity  gradient  is  within  the  range  of  the  commonly  accepted  values.  Finally,  the 
depression  or  shoaling  of  the  interface  results  in  increased  or  decreased  sensitivities, 
respectively,  beneath  the  interface  topography  (Figure  5). 

3.3.  Component-Dependent  Sensitivity  Kernels  and  Utility  of  Three- 
Component  Seismic  Records 

With  the  exception  of  shear-wave  splitting  and  receiver  function  analyses,  the  phase 
or  amplitude  anomaly  of  a  particular  arrival  is  usually  measured  on  only  one  of  the  three- 
component  seismic  records.  Perfectly  good  waveforms  on  the  other  components  are 
often  unutilized.  Using  full  waves  simulated  with  a  finite-difference  method,  we 
demonstrate  that  the  different  components  of  the  same  arrival  at  the  same  receiver  have 
different  traveltime  and  amplitude  sensitivities  to  variations  in  the  velocity  structure 
(Figure  6,  Shen  et  al.,  2008b).  This  is  a  finite-frequency  phenomenon  for  measurements 
derived  from  waveforms.  It  is  important  where  the  scales  of  velocity  heterogeneities  are 
comparable  or  smaller  than  the  width  of  the  Fresnel  zone.  We  calculate  the  Frechet 
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sensitivity  kernels  using  the  scattering-integral  method  (Zhao  et  al.,  2005;  Zhang  et  al., 
2007)  in  conjunction  with  finite-difference  wave  simulation  in  three-dimensional  media. 
The  differences  in  the  sensitivity  kernels  for  the  different  components  vary  with  the  wave 
type,  source-receiver  geometry,  and  source  mechanism  (Shen  et  al.,  2008b).  They  are 
attributed  to  scattered  waves  that  affect  the  waveforms  on  the  different  components  by 
various  amounts.  Thus  the  differential  kernels  between  the  different  components  of  the 
same  arrival  may  enable  us  to  use  the  corresponding  phase  and  amplitude  measurements, 
which  are  relatively  accurate  observations  unaffected  by  uncertainties  in  source  origin 
time  and  location,  to  image  the  earth  structure,  particularly  fine  structures  near  receivers. 

3.4.  Finite-Frequency  Ambient  Noise  Tomography 

Previous  studies  have  demonstrated  the  efficacy  of  extracting  information  on  crustal 
and  shallow  mantle  structure  from  seismic  noise  (e.g.,  Shapiro  and  Campillo,  2004;  Sabra 
et  al.,  2005;  Shapiro  et  al.,  2005;  Yao  et  al.,  2006;  Yang  et  al.,  2007;  Cho  et  al.,  2007;  Lin 
et  al.,  2008).  The  fact  that  the  estimated  Green’s  function  is  based  on  cross-correlation  of 
ambient  noise  makes  it  straightforward  to  adapt  this  type  of  measurements  to  the 
frequency-dependent  phase-delay  anomalies. 

We  have  collected  and  processed  continuous  data  from  50  broadband  stations 
deployed  in  the  Namche  Barwa  seismic  experiment  in  southeast  Tibet  during  2003-2004 
(Sol  et  al.,  2007).  We  use  all  three  components  and  computed  the  cross-correlations 
between  all  possible  pairs  of  components  from  pairs  of  stations  in  order  to  obtain  the 
estimated  Green’s  functions.  The  raw  data  are  cut  into  one-day- length  waveforms.  The 
mean,  trend  and  instrument  response  are  removed  and  the  data  are  bandpass  filtered 
between  4  s  and  150  s  period,  and  one -bit  normalized.  The  data  are  then  cross-correlated 
and  stacked.  The  observed  signals  at  positive  and  negative  and  positive  correlation  lag 
times  correspond  to  waves  propagating  in  opposite  directions  between  two  stations.  We 
measure  the  phase  delays  between  the  estimated  Green’s  function  and  the  synthetics 
calculated  using  a  non-staggered  finite  difference  method  (Zhang  and  Chen,  2006),  which 
takes  into  account  not  only  3D  velocity  heterogeneities  but  also  topography  variations. 
For  the  reference  model,  we  use  a  modified  crustal  model,  CRUST  2.0  to  approximate 
the  3D  wave-speed  variations  in  the  crust,  and  the  akl35  for  the  mantle  part.  We  replace 
the  top  two  sediment  layers  in  CRUST  2.0  with  a  global  sediment  layer  model  digitized 
at  l°xl°  cells.  Topography  is  incorporated  in  the  finite-difference  model. 

Unlike  previous  studies  that  use  the  estimated  Green’s  functions  to  solve  for  phase  or 
group  velocity  maps  and  then  a  point- wise  S  velocity  structure  under  a  ID  assumption 
(e.g.,  Yao  et  al.,  2006;  Yang  et  al.,  2007;  Cho  et  al.,  2007),  we  construct  3D  sensitivity 
kernels  for  the  finite-frequency  Green’s  functions  (Ren  et  al.,  2008).  The  kernels  are 
calculated  using  the  finite-difference  scattering- integral  method  (Zhao  et  al.,  2005;  Zhang 
et  al.,  2007;  Zhang  and  Shen,  2008).  Depending  on  the  wave  period  and  scales  of 
topography,  both  the  shallow  velocity  structure  and  surface  topography  may  significantly 
affect  the  3D  sensitivity  kernels  (Zhang  and  Shen,  2007).  Figure  7  shows  an  example  of 
the  sensitivity  kernels  between  two  stations  in  southeast  Tibet.  The  sensitivity  kernels 
and  the  phase  delays  between  the  estimated  Green’s  functions  and  synthetics  are  inverted 
for  Vp  and  Vs  structure  beneath  southeast  Tibet.  Figure  8  shows  the  Vs  velocity 
perturbations  relative  to  3D  reference  model.  The  broad  low  Vs  anomaly  in  the  lower 
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crust  suggests  that  the  wave  speed  in  the  reference  model  is  too  high  in  the  lower  crust. 
On  the  other  hand,  the  high  Vs  anomaly  in  the  shallow  crust  suggests  that  the  wave  speed 
in  the  reference  model  is  too  low  in  the  shallow  crust. 

4.  DISCUSSION  AND  CONCLUSIONS 

In  the  past  few  years,  the  finite-frequency  theory  has  experienced  a  rapid 
development  from  the  “banana-doughnut”  sensitivity  kernels  based  on  a  ID  reference 
model  and  single-  and  forward-scattering  to  a  full-wave  approach  based  on  a  finite- 
difference  or  spectral-element  method.  Though  approximate,  the  lD-reference  approach 
is  relatively  computationally  efficient,  making  it  well  suited  for  high-frequency 
teleseismic  body  wave  tomography.  The  application  of  this  method  with  finite-frequency 
teleseismic  body  wave  traveltimes  has  led  to  new  observations  of  plume-like  features  in 
the  lower  mantle  (Montelli  et  al.,  2004).  It  has  also  been  adapted  in  various  regional 
crustal  and  mantle  teleseismic  tomography  (Hung  et  al.,  2004;  Yang  et  al.,  2006;  Ren  and 
Shen,  2008;  Liang  et  al.,  2008;  Allen  et  al.,  2008;  Wolfe  et  al.,  2008).  Our  results  show 
that  the  method  works  well  in  places  with  dense  stations  such  as  southeast  Tibet  (Figure 
2).  Although  the  broad  sensitivity  kernels  of  low-frequency  body  waves  help  to  alleviate 
the  problem  of  uneven  sampling  of  the  mantle  structure  beneath  east  Eurasia,  they  do  not 
compensate  for  sparse  station  coverage  and  a  near  vertical  incidence  of  teleseismic  body 
waves  at  shallow  depth  (Figure  3). 

To  date,  two  methods  have  been  developed  to  carry  out  full- wave  tomography 
iteratively  with  3D  reference  models.  One  is  the  adjoint-wavefield  (AW)  method,  which 
back-propagates  the  data  from  the  receivers  to  image  structure  (Tromp  et  al.,  2005;  Liu 
and  Tromp,  2006).  The  other  is  the  scattering-integral  (SI)  method,  which  calculates  and 
stores  the  sensitivity  kernels  for  each  data  functional  (Zhao  et  al.,  2005;  Zhang  et  al., 
2007).  The  AW  approach  has  been  integrated  with  a  spectral-element  method 
(Komatitsch  et  al.,  2005)  for  wave-propagation  simulation  and  the  SI  approach  a  finite- 
difference  method  (Olsen,  1994;  Zhang  and  Chen,  2006).  Each  approach  has  its 
advantages  and  disadvantages.  A  notable  advantage  of  the  SI  method  is  that  it  provides 
both  the  Hessian  and  gradient  of  the  misfit  function  in  the  tomographic  inversion  and 
requires  a  fewer  number  of  simulations  per  model  iteration  to  achieve  the  same  level  of 
data  variance  reduction  (Chen  et  al.,  2007a).  The  disadvantage  is  that  it  requires  a  much 
larger  disk  storage  than  the  AW  method.  Depending  on  the  available  computational 
resources  and  the  numbers  of  sources  and  receivers,  one  approach  may  be  more 
computationally  efficient  than  the  other. 

We  have  focused  our  efforts  on  the  development  and  application  of  the  finite- 
difference,  scattering-integral  method  (FDSIM).  We  have  numerically  validated  the 
scattering-integral  sensitivity  kernels  for  3D  reference  models  (Zhang  et  al.,  2007).  We 
show  that  the  head  wave  arrival-time  and  amplitude  are  influenced  by  the  velocity 
structure  surrounding  the  ray  path  in  a  pattern  that  is  consistent  with  the  Fresnel  zones. 
The  travel-time  sensitivity  of  the  head  wave  is  week  along  the  ray  path  below  the 
interface,  but  non-zero.  Below  the  ray  path,  the  travel-time  sensitivity  decreases  with 
depth  from  the  interface  to  a  minimum  (zero)  before  increasing  to  a  maximum  (absolute 
value),  the  depth  of  which  depends  on  the  wavelength  and  propagation  distance.  We 
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demonstrate  that  the  different  components  of  the  same  arrival  have  different  traveltime 
and  amplitude  sensitivities  to  the  3D  velocity  structure  due  uneven  distribution  of 
scattered  waves  on  the  different  components  of  the  same  arrival  (Shen  et  al.,  2008b). 
Furthermore  we  apply  the  FDSIM  methodology  to  surface-wave  tomography  using 
ambient  seismic  noise.  Unlike  previous  studies  that  use  estimated  Green’s  functions  to 
solve  for  phase  or  group  velocity  maps  and  then  a  point-wise  S  velocity  structure  under  a 
ID  assumption,  we  construct  3D  sensitivity  kernels  for  the  finite-frequency  Green’s 
functions,  taking  into  account  a  3D  reference  model  and  the  effects  of  topographic 
variations.  The  sensitivity  kernels  and  the  phase  delays  between  the  estimated  Green’s 
functions  and  synthetics  are  inverted  for  Vp  and  Vs  structure  beneath  southeast  Tibet 
(Ren  et  al.,  2008).  Results  show  that  the  wave  speed  in  the  reference  model  (a  modified 
CRUST  2.0)  is  too  high  in  the  lower  crust  and  too  low  in  the  shallow  crust  beneath 
southeast  Tibet. 

5.  RECOMMENDATIONS 

The  full- wave  approach  accounts  for  complex  wave  propagation  in  3D  reference 
models,  enables  fuller  utilization  of  an  arrival  on  all  three  components  of  seismic  records, 
and  allows  us  to  linearize  the  inverse  problem  by  iteratively  updating  the  3D  reference 
model.  An  important  benefit  of  physically  realistic  and  accurate  modeling  of  full 
wavefields  in  3D  models  is  the  consistency  of  the  system  of  equations  in  inversion.  This 
is  particularly  important  for  the  integration  of  different  types  of  observations  (P,  S,  and 
surface  waves)  and  physical  properties  (wave  speed  or  elastic  parameters,  attenuation, 
and  anisotropy)  in  inversion,  which  is  a  challenging  but  essential  step  to  obtain  a  coherent 
and  self-consistent  image  of  the  crust  and  mantle.  One  notable  benefit  of  the  full-wave 
approach  is  that  the  new  method  makes  it  unnecessary  to  identify  and  separate  body  and 
surface  waves.  It  thus  provides  a  natural  way  to  integrate  the  constraints  from  various 
wave  types.  With  the  ever-increasing  high-performance  computational  power,  this  full- 
wave  approach  will  become  a  more  and  more  practical  and  powerful  tool.  Peta-scale 
computing  systems  (e.g.,  Blue  Waters,  a  NSF-sponsored  system),  for  instance,  will  be 
available  for  scientific  applications  in  the  next  few  years.  We  anticipate  that  such 
facilities  will  take  the  full- wave  approach  to  a  new  level. 

The  applications  of  our  FDSIM  are  limited  by  computational  resources.  The  greater 
the  computational  resources  we  have,  the  larger  the  geographic  area  we  can  study  or  the 
higher  frequency  waves  we  can  use.  Nevertheless,  to  utilize  high-frequency  (0.2  -  1  Hz) 
teleseismic  body  waves  in  regional  tomography,  it  will  be  necessary  in  the  near  future  to 
develop  a  coupled  1D-3D  approach,  in  which  waves  outside  of  the  volume  of  interest 
propagate  within  a  ID  model  and  those  inside  the  volume  of  interest  are  treated  with  a 
full  3D  numerical  method. 
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Figure  1  Topography  of  Eastern  Eurasia  and  the  distribution  of  broadband  seismic 
stations  used  in  this  study  (triangles).  The  box  outlines  a  detailed  finite-frequency 
teleseismic  traveltime  tomographic  study  in  southeastern  Tibet. 
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Figure  2  (a)  P-velocity  perturbations  at  several  depths  beneath  southeastern  Tibet:  75  km, 
1 12  km,  150  km,  262  km,  375  km,  and  450  km.  Black  lines  in  the  75-km  horizontal  slice 
represent  the  known  faults.  Notice  the  north-south  trending  rift  coincides  with  a 
prominent  north-south  trending  low- velocity  anomaly.  Line  AB  on  the  1 12-km-depth 
image  marks  the  location  of  the  vertical  profile  in  (c).  (b)  S-velocity  perturbations  at  the 
same  depths,  (c)  Cross-section  A-B  through  P  and  S-wave  tomographic  models,  along 
the  Indus-Yarlung  suture  and  across  a  rift  in  southeastern  Tibet.  From  Ren  and  Shen 
(2008). 
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Figure  3  The  P-velocity  variations  relative  to  iasp91  at  180,  270  and  540  km  depth. 
Triangles  in  the  upper  left  panel  show  the  distribution  of  the  broadband  stations  used  in 
eastern  Eurasia. 
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Figure  4  (a)  The  traveltime  sensitivity  kernels  for  the  head  wave  recorded  on  the  vertical 
component  on  the  horizontal  planes  at  depths  below  (21.6  km,  24  km,  26.4  km,  28.8  km) 
and  above  (19.2  km)  the  interface  that  separates  the  low  and  high  velocity  layers.  The 
green  line  marks  the  head  wave  ray  path.  The  source  is  located  on  the  left  side  of  this 
figure.  The  negative  (red  colors)  and  positive  (blue  colors)  values  are  so  defined  that  a 
low-velocity  anomaly  located  in  the  region  of  the  negative  kernels  results  in  a  travel  time 
delay  and  the  same  velocity  perturbation  in  the  region  of  positive  kernels  leads  to  an 
earlier  head  wave  arrival,  (b)  The  traveltime  sensitivity  kernel  for  the  head  wave  on  a 
vertical  profile  half  way  between  the  source  and  receiver  and  perpendicular  to  the  ray 
path.  The  dashed  line  is  the  layer  interface,  (c)  The  traveltime  sensitivity  kernels  for  the 
head  wave  on  the  vertical  profile  containing  the  source  and  receiver.  The  full  range  of  the 
kernels  is  -4.83-1.655xl0"n  s/m3.  From  Zhang  et  al.  (2007). 
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Figure  5  The  kernels  for  the  vertical  component  of  head  waves  in  the  models  with  up- 
bending  and  down-bending  layer  interface  between  the  source  and  receiver.  Except  for 
the  perturbed  interface,  all  other  parameters  of  the  models  are  the  same  as  those  of  the 
model  in  Fig.  4.  Note  that  the  interface  variation  is  cylindrical  with  the  symmetric  axis 
perpendicular  to  the  plane  containing  the  source-receiver  path,  not  spherical.  The  dashed 
line  indicates  the  layer  interface  and  the  green  line  highlights  the  head  wave  ray  path.  The 
meanings  of  the  positive  and  negative  values  are  the  same  as  those  in  Fig.  4.  (a)  A 
cylindrical  down-bending  interface  with  the  maximum  depression  of  10  km  and  a 
horizontal  dimension  of  35  km.  The  full  range  of  the  kernel  is  -5.19-1.75x10'  's/m3  (b) 
An  up-bending  interface.  The  maximum  shoaling  is  5  km,  and  the  horizontal  dimension 
of  the  interface  perturbation  is  35  km.  The  full  range  of  the  kernel  is  -3.82-1.43x10" 
ns/m3.  From  Zhang  et  al.,  (2007). 
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Figure  6  Delay-time  (left)  and  amplitude  (right)  sensitivity  kernels  for  the  z  and  x 
components  and  their  differential  kernel.  The  surface  receiver  (30.8  km,  0)  is  at  a 
relatively  short  horizontal  distance  from  the  explosive  source  at  depth  (22.8  km,  30.0 
km).  The  negative  (red  colors)  and  positive  (blue  colors)  values  for  traveltimes  are  so 
defined  that  a  low-velocity  anomaly  located  in  the  region  of  the  negative  kernels  results 
in  a  phase  delay  and  the  same  velocity  perturbation  in  the  region  of  positive  kernels  leads 
to  an  earlier  arrival  as  measured  by  waveform  cross-correlation.  For  amplitude,  the 
negative  (red  colors)  and  positive  (blue  colors)  values  are  so  defined  that  a  high-velocity 
anomaly  located  in  the  region  of  the  negative  kernels  results  in  an  amplitude  reduction 
and  the  same  velocity  perturbation  in  the  region  of  positive  kernels  leads  to  an  increase  in 
amplitude.  From  Shen  et  al.  (2008). 
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Figure  7  Cross-sections  of  phase-delay  sensitivities  (in  s/m3)  to  Vs  perturbations  between 
the  vertical  components  of  two  seismic  stations  in  southeast  Tibet  for  wave  periods  of  20- 
30  s.  The  phase-delays  are  sensitive  also  to  Vp  perturbations  at  shallow  depths  (not 
shown). 
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Figure  8  Preliminary  Vs  perturbation  models  from  the  inversion  of  estimated  Green’s 
function  phase  delays  at  12  km  (upper  panel)  and  52  km  depth  (lower  panel).  The 
anomaly  is  relative  to  the  3D  reference  model,  a  modified  CRUST  2.0.  The  color  scale  is 
±15%  (top  panel)  and  ±5%  (lower  panel),  respectively. 
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AFRL 

Air  Force  Research  Laboratory 

AW 

Adjoint  wavefield 

DMC 

Data  Management  Center 

EHB 

Engdahl-van  der  Hilst-Buland 

FDSIM 

finite-difference,  scattering-integral  method 

GSN 

Global  Seismographic  Network 

IRIS 

Incorporated  Research  Institutions  for  Seismology 

JISNET 

Japan-Indonesia  Seismic  Network 

NEIC 

National  Earthquake  Information  Center 

PASSCAL 

Program  for  Array  Seismic  Studies  of  the  Continental  Lithosphere 

SI 

Scattering  Integral 
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Short  Note 

Frequency-Dependent  Crustal  Correction  for  Finite-Frequency 

Seismic  Tomography 

by  Ting  Yang*  and  Yang  Shen 


Abstract  Removing  the  crustal  signature  from  teleseismic  travel  times  is  an  im¬ 
portant  procedure  to  reduce  the  trade-off  between  crustal  and  mantle  velocity  het¬ 
erogeneities  in  seismic  tomography.  Because  reverberations  of  long-  and  short-period 
body-wave  arrivals  in  the  crust  affect  the  waveforms  of  the  direct  arrivals  differently, 
the  crustal  effects  on  travel  times  measured  by  waveform  cross  correlation  are  fre¬ 
quency  dependent.  With  synthetic  responses  of  selected  crustal  models,  this  short 
note  illustrates  the  significance  of  frequency-dependent  crustal  corrections  to  finite- 
frequency  body-wave  travel-time  tomography.  The  differences  in  crustal  correction 
between  long-  and  short-period  body  waves  at  the  same  station  can  be  as  large  as 
0.6  sec,  depending  on  the  crustal  thickness,  velocity  contrast  at  the  Moho,  and  lay¬ 
ering  within  the  crust. 


Introduction 

Variations  in  station  elevation  and  cmstal  structure  be¬ 
neath  seismic  stations  are  responsible  for  a  significant  frac¬ 
tion  of  observed  teleseismic  body-wave  travel-time  anoma¬ 
lies.  In  regions  having  large  variations  in  crustal  thickness 
(e.g.,  eastern  Eurasia  and  central  Andes)  or  complex  cmstal 
structure,  the  crust  may  account  for  half  of  the  observed 
travel-time  residuals  (Kissling,  1993;  Beck  etal. ,  1996;  Mar¬ 
tin  et  al ,  2005).  In  teleseismic  body- wave  tomographic  stud¬ 
ies,  the  cmstal  and  shallow-mantle  structure  usually  cannot 
be  directly  resolved  in  the  inversion  because  incident  waves 
from  teleseismic  events  propagate  nearly  vertically  and  do 
not  cross  each  other  at  shallow  depths  (e.g.,  Wolfe  et  al ., 
1997).  If  cmstal  anomalies  are  not  removed  from  travel-time 
residuals,  they  tend  to  be  mapped  into  deeper  mantle,  caus¬ 
ing  artificial  velocity  anomalies  in  tomographic  models  (e.g., 
Waldhauser  et  al ,  2002).  Travel-time  correction  for  shallow 
structure,  therefore,  has  been  an  essential  procedure  in  seis¬ 
mic  tomography  to  improve  the  resolution  in  the  mantle. 

There  are  two  common  approaches  to  cope  with  this 
problem.  The  first  is  to  add  an  additional  unknown  for  each 
station  to  be  solved  in  the  inversion  (e.g.,  Evans  and 
Achauer,  1993;  VanDecar  et  al ,  1995;  Wolfe  et  al ,  1997; 
Foulger  et  al ,  2001).  These  additional  unknowns,  or  station 
terms,  are  used  to  account  for  travel-time  anomalies  asso¬ 
ciated  with  the  cmstal  and  shallow-mantle  structure  and  are 
solved  simultaneously  with  the  velocity  structure  in  the  in¬ 
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version.  One  notable  characteristic  of  this  approach  is  that 
the  station  terms  include  the  effects  from  both  the  crust  and 
shallow  mantle.  Thus  for  studies  with  a  large  station  spacing, 
the  shallow-mantle  structure  is  usually  “absorbed”  into  the 
station  terms  and  unresolved.  By  definition,  the  station  terms 
do  not  vary  with  frequencies,  epicentral  distances,  and  back- 
azimuths  of  incoming  waves. 

For  regions  where  the  cmstal  structure  has  been  deter¬ 
mined  from  other  observations  such  as  surface- wave  studies, 
receiver  function  analyses,  and/or  reflection  and  refraction 
seismic  profiles,  the  cmstal  effects  relative  to  a  reference 
model  can  be  removed  by  subtracting  the  travel-time  anom¬ 
alies  accrued  within  the  crust  from  the  total  travel-time  de¬ 
lays.  The  residual  travel-time  delays  are  then  used  to  invert 
the  mantle  structure  (e.g.,  Dawson  et  al ,  1990;  Allen  et  al , 
2002;  Keyser  et  al ,  2002).  Compared  with  the  station-term 
method,  this  approach  is  valuable  for  imaging  the  shallow- 
mantle  structure  (e.g.,  Allen  et  al ,  2002),  in  particular,  when 
the  a  priori  crustal  structure  is  considered  to  be  well  con¬ 
strained.  Crustal  corrections  in  this  approach  may  vary  with 
epicentral  distances  and  backazimuths  of  incoming  waves. 
They  are  usually  calculated  for  the  direct  arrival  under  ray 
theory  without  any  consideration  to  the  reverberation  of  the 
wave  in  the  crust  (e.g.,  Allen  et  al ,  2002;  Waldhauser  et  al , 
2002;  Hung  et  al ,  2004;  Martin  et  al ,  2005).  When  the 
arrival  time  of  a  seismic  phase  is  handpicked  from  the  onset 
of  the  arrival,  ray  theory  is  commonly  used  and  justified  by 
the  argument  that  the  onset  represents  the  arrival  of  the 
highest-frequency  signal  in  the  seismogram.  But  it  is  often 
difficult  to  pick  the  onset  of  the  arrival  accurately  when  the 
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arrivals  are  long  period,  emergent,  or  of  less-than-excellent 
signal-to-noise  ratio. 

The  availability  of  digital  seismic  records  has  made  it  a 
common  practice  to  obtain  travel  times  by  cross-correlating 
seismic  waveforms  (VanDecar  and  Crosson,  1990).  For 
finite-frequency  seismic  signals,  reverberations  of  the  arrival 
in  the  crust  beneath  a  station  may  overlap  with  the  direct 
arrival  and  thus  may  significantly  change  the  waveform  of 
the  direct  arrival.  The  travel  times  obtained  by  waveform 
cross-correlation  are  therefore  sensitive  to  crust  reverbera¬ 
tions,  especially  at  low  frequencies  and  for  the  crust  with 
strong  reflectors  at  shallow  depth.  As  a  result,  the  correction 
calculated  for  the  direct  ray  path  under  ray  theory  cannot 
fully  account  for  the  effects  of  crustal  variations  on  travel 
times  measured  by  waveform  cross  correlation. 

This  problem  is  often  neglected  in  ray-based  tomogra¬ 
phy  and  has  now  become  more  severe  for  the  newly  devel¬ 
oped  finite-frequency  seismic  tomography  (FFST),  in  which 
the  travel-time  calculation  is  based  on  3D  sensitivity  kernels 
(Dahlen  et  al. ,  2000)  instead  of  ray  theory.  Because  the  3D 
volumes  of  low-frequency  kernels  are  broader  than  those  of 
higher-frequency  signals,  low-frequency  signals  offer  con¬ 
straints  on  the  structure  farther  away  from  geometrical  rays. 
This  more  realistic  representation  of  wave  propagation  helps 
to  obtain  a  smoother,  more  even  sampling  of  the  velocity 
structure  beneath  stations.  At  ocean  bottom  and  ocean  island 
environments,  where  the  microseism  is  high  (Wilcock  et  al ., 
1999),  most  useful  arrivals  with  good  signal-to-noise  ratios 
are  in  the  low-noise  notch  (roughly  0.03-0.1  Hz;  Webb, 
1998),  the  low-frequency  range  in  teleseismic  body-wave 
tomography.  To  isolate  the  microseism,  Hung  et  al.  (2004), 
Shen  and  Hung  (2004),  and  Yang  etal.  (2006)  filtered  broad¬ 
band  P-wave  records  in  the  high-,  intermediate-,  and  low- 
frequency  bands  (0.5-2,  0. 1-0.5,  and  0.03-0.1  Hz,  respec¬ 
tively),  and  combined  the  travel  times  measured  in  the  three 
frequency  ranges  in  joint  tomographic  inversions.  About 
60%  of  the  useful  P-wave  arrivals  in  the  Iceland  data  set, 
for  example,  is  long  period.  For  these  reasons,  long-period 
travel-time  data  and  their  correction  for  crustal  effects  are 
important  in  FFST,  which  offers  the  potential  to  utilize  more 
information  in  broadband  waveforms  than  ray  theory. 

Here  we  present  a  crustal  correction  method,  in  which 
crustal  reverberations  for  arrivals  with  different  frequencies 
are  taken  into  account.  Using  synthetic  waveforms,  we  show 
that  the  crustal  correction  can  be  approximated  by  the  wave¬ 
form  shift  measured  by  cross  correlation  of  the  response  of 
a  given  crustal  model  and  that  of  the  reference  model  in  a 
narrow-frequency  band. 

Effects  of  Crustal  Reverberations  on  Travel  Times 

We  use  synthetic  seismograms  to  demonstrate  that  the 
effects  of  body-wave  reverberations  in  the  crust  beneath  a 
station  on  travel  times  measured  by  waveform  cross  correla¬ 
tion  are  frequency  dependent.  We  construct  synthetic  seis¬ 
mograms  by  convolving  a  source  time  series  with  the  impulse 


response  of  a  crustal  model.  Our  source  function  is  chosen  to 
be  the  first  derivative  of  a  Gaussian  pulse,  of  the  form 

src(t)  =  -8 7i2  ft  -  i  rj 


where  the  frequency  content  of  this  source  function  is  con¬ 
trolled  by  r.  This  time  function  is  a  full-cycle  pulse  (shown 
in  Fig.  la).  The  impulse  responses  are  computed  for  a  given 
layered  crustal  model  and  a  ray  parameter  using  Kennett’ s 
reflectivity  matrix  approach  (Kennett,  1983).  The  calculation 
yields  a  full  response  including  mode  conversions.  For  our 
purposes,  we  consider  the  response  at  the  free  surface. 

Figure  1  shows  the  procedure  of  calculating  synthetic 
seismograms  from  two  crustal  models  and  the  frequency- 
dependent  crustal  corrections  obtained  from  waveform  cross 
correlation.  One  of  the  models  is  the  IASP91  (Kennett  and 
Engdahl,  1991).  The  other  is  the  crustal  structure  beneath 
the  Global  Seismic  Network  (GSN)  station  DGAR  (7.412°  S, 
72.453°  E),  an  ocean  island  site  with  a  relatively  thin  crust, 
from  CRUST2.0  (Bassin  et  al .,  2000).  For  simplicity,  the 
first  four  layers  in  CRUST2.0  (water,  ice,  and  soft  and  hard 
sediment)  are  not  included  in  the  calculations.  This  is  rea¬ 
sonable  because  permanent  seismic  stations  are  usually  set 
up  on  bedrocks.  The  thickness,  Vp,  Vs,  and  density  of  each 
layer  in  the  two  crustal  models  are  shown  in  Figure  lb.  Fig¬ 
ure  lc  shows  the  calculated  impulse  responses  of  these  two 
models  to  an  incident  P  wave  by  the  reflectivity  matrix 
method  (Kennett,  1983).  The  impulse  responses  are  then 
convolved  with  the  source  time  function  (Fig.  la)  to  generate 
synthetic  seismograms  (Fig.  Id).  To  find  the  frequency  de¬ 
pendence  of  the  travel-time  differences  between  the  two 
models,  we  filtered  these  seismograms  using  the  high-, 
intermediate-,  and  low-frequency  bands  as  in  Hung  et  al. 
(2004)  (0.5-2,  0. 1-0.5,  and  0.03-0.1  Hz,  respectively). 
Relative  travel-time  delays  or  crustal  corrections  for  the 
three  frequency  bands  are  then  measured  from  the  cross 
correlation  of  the  two  filtered  seismograms.  As  shown  in 
Figure  le,  while  the  travel-time  differences  are  nearly  iden¬ 
tical  for  high-  and  intermediate-frequency  signals,  the  arrival 
at  DGAR  at  the  low  frequency  is  much  earlier  (by  an  addi¬ 
tional  0.65  sec)  than  at  the  intermediate  and  high  frequen¬ 
cies,  due  to  the  overlap  of  crustal  reverberations  with  the 
long-period  direct  arrivals  and  the  smaller  separation  be¬ 
tween  the  reverberations  and  the  direct  arrival  at  the  station 
above  a  thin  crust. 

To  see  how  variations  in  crustal  thickness  alone  affect 
the  travel  times  measured  by  waveform  cross  correlation,  we 
systemically  generated  crustal  models  with  different  crustal 
thicknesses  and  repeated  the  preceding  procedure  to  calcu¬ 
late  the  corrections  relative  to  the  IASP91  model.  In  this 
experiment,  the  crust  is  a  two-layer  model,  in  which  the 
thickness  of  each  layer  is  increased  or  decreased  by  the  same 
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Figure  1 .  Procedure  to  generate  synthetic  seismograms  on  the  vertical  component  for  frequency-dependent 
crustal  correction,  (a)  Source  time  function.  The  power  spectrum  of  this  example  peaks  at  about  1  Hz  and  has 
energy  between  0.001  and  2  Hz.  (b)  The  IASP91  crustal  model  and  a  three-layer  crustal  model  beneath  the 
DGAR  station  from  CRUST2.0  (Bassin  et  al. ,  2000).  Values  within  layers  are  thickness,  Vp,  Vs,  and  density, 
respectively,  (c)  Responses  of  the  IASP91  model  and  the  DGAR  crustal  model.  The  major  reverberation  phases 
are  marked.  PdP  ( PdS )  stands  for  a  P  wave  reflected  at  the  surface  and  then  reflected  (converted  to  S)  at  the 
discontinuity  d  (1,  2,  or  the  Moho).  (d)  Synthetic  seismograms  obtained  from  the  convolution  of  the  source 
time  function  and  the  crustal  responses,  (e)  Bandpass-filtered  seismograms  for  the  IASP91  (dashed)  and  DGAR 
crust,  and  travel-time  shifts  between  them  measured  by  cross  correlations.  (« continued ) 
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Figure  1.  (continued),  “f”  marks  the  peak  of  seismograms  from  the  IASP91  model,  which  is  also  the 
cross-correlation  reference  time;  “cc”  is  the  travel-time  shift  of  the  two  seismograms  measured  from  cross 
correlation  (i.e.,  crustal  correction).  Shaded  portions  with  dashed  (IASP91)  and  solid  (DGAR)  frames 
indicate  the  windows  for  cross  correlation.  The  crustal  correction  from  the  ray  theory  is  0.64  sec,  and  the 
ray  parameter  in  this  experiment  is  0.06  sec/km.  The  difference  between  the  value  predicted  by  the  ray 
theory  and  those  obtained  from  waveform  crosscorrelation  is  attributable  to  the  large  reverberation  from 
the  shallow  crust  beneath  DGAR  {PIP). 


amount  incrementally  while  all  other  parameters  are  kept  the 
same  as  those  in  IASP91.  Figure  2  shows  the  relationship 
between  the  crustal  thickness  and  corrections.  Corrections 
for  low-,  intermediate-,  and  high-frequency  signals  differ 
when  the  crustal  thickness  is  less  than  that  of  IASP91 
(35  km),  with  the  most  significant  differences  in  the  thick¬ 
ness  range  of  6-24  km.  For  the  crust  with  a  thickness  near 
or  greater  than  35  km,  the  separation  between  the  main  re¬ 
verberation  phase  {PmP)  and  the  direct  arrival  is  greater  than 
1 1  sec,  sufficiently  large  that  the  reverberation  does  not  over¬ 
lap  with  the  direct  arrival  and  thus  does  not  affect  the  mea¬ 
sured  travel  times. 

To  see  how  internal  layering  affects  the  travel  time  mea¬ 
sured  by  waveform  cross  correlation,  we  did  a  similar  ex¬ 
periment  with  a  more  complex,  three-layer  crustal  model 
(Fig.  3).  Velocities  and  densities  for  this  crust  model  are 
derived  from  the  average  values  of  the  crustal  structure  be¬ 
neath  180  eastern  Eurasia  stations  from  CRUST2.0  (water, 
ice,  and  sediment  layers  are  also  excluded),  though  we  ar¬ 
bitrarily  fix  the  thickness  of  each  layer  to  be  one  third  of  the 
total  crustal  thickness.  These  parameters  are  shown  in 
Table  1.  Unlike  in  Figure  2,  the  frequency  dependence  in 
travel  times  does  not  disappear  when  the  crustal  thickness  is 
equal  to,  or  larger  than  that  of  the  reference  model  (35  km). 
This  is  attributable  to  reverberations  in  the  shallow  crust  in 
the  three-layer  model.  So  in  addition  to  crustal  thickness,  the 


internal  layering  in  the  crust,  in  particular,  the  strong  reflec¬ 
tion  in  the  shallow  crust,  may  also  have  a  significant  effect 
on  the  waveforms  and  thus  the  travel  times  measured  by 
waveform  cross  correlation. 

The  preceding  synthetic  experiments  illustrate  that,  for 
travel  times  measured  by  waveform  cross  correlation,  the 
effects  of  crustal  reverberations  are  frequency  dependent  and 
must  be  taken  into  consideration  in  tomographic  inversions. 
Crustal  corrections  calculated  simply  from  ray  theory  for  the 
direct  arrival  may  leave  large  biases  in  the  travel  times  of 
low-frequency  waves  and  even  short-period  waves  if  strong 
reflectors  exist  in  the  shallow  (depth  <5  km)  crust.  Since 
most  of  the  oceanic  crust  and  a  significant  portion  of  the 
continental  crust  have  a  thickness  within  the  range  of  5  to 
30  km,  this  bias  in  travel  times  may  significantly  degrade 
the  resolution  of  the  mantle  structure  in  seismic  tomography 
(van  der  Hilst  and  de  Hoop,  2005).  Frequency-dependent 
crustal  corrections  should  be  applied  in  FFST  and  other  stud¬ 
ies  using  low-frequency  seismic  records. 

Practical  Considerations 

If  an  accurate  crustal  structure  is  available,  the  effect  of 
crustal  reverberations  on  waveforms  can  be  removed  by  de¬ 
convolving  the  crustal  response  from  the  observed  seismo¬ 
grams.  The  crustal  corrections  calculated  from  ray  theory 
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Figure  2.  Effect  of  the  crustal  thickness  on  crustal 
corrections  for  the  high-  (0.5-2.0  Hz),  intermediate- 
(0. 1-0.5  Hz),  and  low-  (0.03-0.1  Hz)  frequency  sig¬ 
nals.  The  crust  is  a  two-layer  model  with  the  same 
Vp,  Vs,  and  density  as  those  of  IASP91.  The  thick¬ 
nesses  of  the  upper  and  lower  crust  are  increased  or 
decreased  by  the  same  amount  in  each  calculation. 
The  source  time  function  is  shown  in  Figure  la. 


can  then  be  used  properly  for  signals  at  all  frequencies.  How¬ 
ever,  this  approach  may  not  be  practical  since  the  decon¬ 
volution  process  is  sensitive  to  the  response,  which  requires 
a  very  accurate  crustal  model  (Obayashi  et  al. ,  2004). 

Crustal  corrections  can  also  be  approximated  without 
the  deconvolution  of  recorded  waveforms  by  following  the 
steps  in  Figure  1 .  In  practice,  this  may  also  be  achieved  by 
calculating  the  impulse  response  of  the  crustal  model  di¬ 
rectly,  and  filtering  the  response  within  the  specified  fre¬ 
quency  bands.  Real  seismic  records  may  have  a  higher  en¬ 
ergy  toward  the  low-frequency  end  of  the  selected  frequency 
band  or  vice  versa.  Thus  real  data  may  have  a  spectrum  that 
is  somewhat  different  from  the  impulse  response,  whose 
power  spectrum  is  flat.  This  mismatch  between  the  spectra 
of  real  records  and  the  impulse  response  may  introduce  er¬ 
rors  into  the  crustal  corrections  for  travel  times  measured  by 
waveform  cross  correlation.  Our  synthetic  experiments 
show,  however,  that  this  error  is  likely  negligible  with  nar¬ 
row  frequency  bands.  Figure  4  shows  the  difference  between 
the  correction  obtained  from  synthetic  seismograms  with  a 
known  source  time  function  (ccl)  and  that  directly  from  the 
impulse  response  of  the  crust  (cc2).  We  conducted  this  syn¬ 
thetic  experiment  for  two  sets  of  crustal  models  and  several 
different  source  time  functions  having  various  power  spec¬ 
tra.  For  the  three  frequency  bands  and  all  crustal  models 
with  different  thicknesses,  most  of  the  differences  between 
the  two  approaches  are  less  than  one  sampling  interval 


Figure  3.  Effects  of  the  crustal  thickness  and  the 
internal  layering  on  crustal  corrections  for  the  high- 
(0.5-2.0  Hz),  intermediate-  (0. 1-0.5  Hz),  and  low- 
(0.03-0.1  Hz)  frequency  signals.  The  crust  is  a  three- 
layer  model  with  Vp,  Vs,  and  density  being  the  mean 
values  of  the  crust  beneath  seismic  stations  in  eastern 
Eurasia  from  CRUST2.0  (layer  5,  6,  and  7)  (Bassin 
et  al .,  2000).  The  thickness  of  each  layer  is  one  third 
of  the  total  crustal  thickness.  The  average  P-wave 
crustal  velocity  (6.55  km/sec)  is  higher  than  that  in 
IASP91  (6.08  km/sec).  The  source  time  function  is 
shown  in  Figure  la. 


Table  1 

Average  Crustal  Structure  beneath  Seismic  Stations  in  Eastern 
Eurasia*  from  CRUST2.0 


Vp  (km/sec) 

Vs  (km/sec) 

Density  (g/cm3) 

Upper  crust 

5.97 

3.39 

2.71 

Middle  crust 

6.55 

3.67 

2.88 

Lower  crust 

7.13 

3.95 

3.06 

*  Stations  are  within  the  range  of  (20°  S  to  60°  S,  60°  E  to  160°  E). 


(0.05  sec),  and  the  biggest  one  is  0.08  sec.  This  result  sug¬ 
gests  that,  to  first  order,  the  crustal  corrections  within  indi¬ 
vidual  frequency  bands  are  not  very  sensitive  to  the  spectrum 
of  the  incoming  wave  and  one  may  choose  to  use  corrections 
measured  directly  from  impulse  responses  to  approximate 
those  for  observed  waveforms. 

With  the  preceding  approximations,  a  first-order  crustal 
correction  can  be  obtained  through  the  following  steps: 

1.  Calculate  the  ray  parameter  based  on  the  station  and 
earthquake  locations. 

2.  Calculate  the  responses  of  the  crustal  model  beneath  the 
station  and  the  reference  model  for  the  ray  parameter. 
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Figure  4.  Differences  between  crustal  corrections  from  synthetic  seismograms  (cel) 
and  those  from  impulse  responses  (cc2)  for  the  high-  (0.5-2.0  Hz),  intermediate-  (0.1- 
0.5  Hz),  and  low-  (0.03-0.1  Hz)  frequency  bands,  (a)  The  two-layer,  IASP91-type  crust, 
(b)  The  three-layer  crust,  as  shown  in  Table  1.  “T”  is  the  apparent  period  of  the  source 
time  source  (t)  in  equation  (1). 


3.  Filter  the  responses  in  the  same  frequency  ranges  as  those 
of  real  data. 

4.  Cross -correlate  the  two  filtered  responses  from  the  ref¬ 
erence  model  and  the  model  for  a  station  to  measure  the 
correction  for  a  specified  frequency  band. 

5.  Apply  the  correction  to  the  corresponding  travel  time  of 
the  observed  waveform  at  the  station. 

To  get  a  quantitative  estimate  of  the  magnitude  of 
frequency-dependent  crustal  corrections  in  a  real-world  sit¬ 
uation,  we  applied  this  crust-correction  procedure  to  broad¬ 
band  seismograms  recorded  at  stations  in  the  western  Pacific. 
The  crustal  models  beneath  these  stations  are  from 
CRUST2.0  (Bassin  et  al. ,  2000).  Figure  5  shows  an  example 
of  the  effect  of  the  frequency-dependent  crustal  correction 
on  travel  time  shifts,  in  which  the  magnitude  of  crustal  cor¬ 
rection  is  comparable  to  the  relative  delays  measured  by 
cross  correlation  and  the  difference  between  the  high-  and 
low-frequency  corrections  is  as  much  as  25%  of  the  travel¬ 
time  delays. 

There  are  two  other  assumptions  in  the  previous  ex¬ 
amples  that  one  may  consider  in  the  crustal  correction.  First, 
the  uppermost  mantle  velocity  affects  the  impedance  con¬ 
trast  at  the  Moho  and  thus  has  an  effect  on  the  waveform  of 
the  crustal  reverberation  and  on  the  crustal  correction.  In 


many  cases  (e.g.,  the  global  crustal  model  CRUST2.0),  the 
uppermost  mantle  velocity  is  unavailable.  One  has  to  assume 
a  velocity  for  the  uppermost  mantle.  To  assess  how  this  as¬ 
sumption  affects  the  crustal  correction,  we  changed  the  up¬ 
permost  mantle  velocity  and  compared  crustal  corrections 
with  that  based  on  the  IASP91  velocity.  A  ±3.7%  variation 
in  the  uppermost  mantle  velocity  (7.75-8.35  km/sec)  causes 
a  very  small  difference  (0.035  sec)  in  the  measured  low- 
frequency  travel  times.  Therefore,  to  first  order,  the  effect  of 
the  uppermost  mantle  velocity  is  negligible. 

Second,  we  assume  flat  crustal  interfaces  in  these  cal¬ 
culations  of  the  crustal  response.  The  real  crust  has  2D  or 
3D  variations  in  velocity  and  interfaces,  which  have  more 
complex  influences  on  the  waveforms  of  recorded  broad¬ 
band  signals.  For  regions  having  well-constrained  2D  or  3D 
crustal  models,  the  crustal  responses  can  be  calculated  with 
other  techniques  (e.g.,  Komatitsch  and  Vilotte,  1998;  Fred- 
eriksen  and  Bostock,  2000;  Kang  and  Baag,  2004). 

This  paper  focuses  on  P  waves.  We  note  that  teleseismic 
S  waves  are  usually  observed  at  periods  of  several  seconds 
and  longer.  Depending  on  the  frequency  of  the  S  arrivals  and 
the  crustal  structure,  reverberations  of  S  waves  in  the  crust 
may  also  overlap  with  the  direct  S  arrival  and  thus  affect  the 
travel  times  measured  by  waveform  cross  correlation.  The 
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Figure  5.  The  filtered  high-  (0.5-2.0  Hz)  (a)  and  low-  (0.03-0.1  Hz)  (b)  frequency 
seismograms  at  six  stations  in  the  western  Pacific.  Seismograms  are  aligned  according 
to  the  time  shifts  determined  from  waveform  cross  correlation.  The  vertical  bars  mark 
the  time  shifts  (relative  to  the  “r”  bar  in  the  first  trace)  calculated  from  the  waveform 
cross  correlation  (“c”).  “h”  and  “1”  represent  the  adjusted  time  shifts  after  the  crustal 
corrections  estimated  from  the  impulse  response  in  the  high-  and  low-frequency  bands, 
respectively.  In  low-frequency  seismograms  (b),  the  time  shifts  with  crustal  corrections 
calculated  from  the  high-frequency  impulse  response  (“h”)  are  also  shown  for  com¬ 
parison  with  “1”.  Values  after  station  names  are  the  thickness  of  the  crust  beneath  the 
station  in  kilometers  (Bassin  et  ah ,  2000)  and  the  ray  parameter  in  seconds  per  kilometer 
The  earthquake  occurred  at  02:45:59  on  7  September  2001,  and  the  epicenter  was  at 
(13.16°  S,  97.30°  E). 


frequency-dependent  S  crustal  correction  can  be  determined 
in  the  same  way  as  outlined  for  P  waves. 

Conclusions 

Through  a  series  of  synthetic  experiments,  we  demon¬ 
strate  that  the  effects  of  crustal  structure  on  travel  times  mea¬ 
sured  by  waveform  cross  correlation  are  frequency  depen¬ 
dent.  The  difference  in  crustal  correction  between  high-  and 
low-frequency  arrivals  could  reach  as  large  as  0.65  sec  for 
a  relatively  thin  crust  and  crust  with  strong  shallow  rever¬ 
berations.  As  a  result,  this  frequency  dependence  in  travel 
times  measured  by  waveform  cross  correlation  cannot  be 
ignored  in  finite-frequency  tomography  or  any  tomography 
using  low  frequency  or  broadband  body  waves.  The 
frequency-dependent  crustal  correction  can  be  approxi¬ 
mated,  to  first  order,  by  cross-correlating  the  impulse  re¬ 
sponses  of  a  crust  model  filtered  in  a  narrow  frequency  band. 
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[i]  While  several  mechanisms  have  been  suggested  to  explain  the  evolution  of  the 
Tibetan  Plateau,  observational  constraints  on  the  deep  lithospheric  processes  have  been 
sparse,  and  previous  seismic  studies  were  mostly  along  profiles  perpendicular  to  the 
collision  front  of  the  Indian  and  Eurasian  plates.  In  this  study,  we  show  tomographic 
evidence  for  the  delamination  of  the  mantle  lithosphere  beneath  southeastern  Tibet,  a 
process  in  which  the  entire  mantle  lithosphere  peels  away  from  the  crust  along  the  Moho 
and  thus  is  a  mechanism  for  rapid  thinning  of  the  lithosphere.  Our  P  and  S  wave  velocity 
models  show  the  presence  of  a  low-velocity  anomaly  in  the  crust  and  upper  mantle 
down  to  ~300  km  depth  beneath  a  north- south  trending  rift  zone  in  southeastern  Tibet. 

This  low-velocity  anomaly  is  situated  above  a  tabular,  high-dipping-angle,  high-velocity 
anomaly  that  extends  into  the  upper  mantle  transition  zone.  The  VP/VS  ratio  of  this 
high-velocity  anomaly  suggests  that  temperature  variations  are  not  the  only  cause  of  the 
anomaly  and  a  highly  melt-depleted  mantle  is  required.  These  observations  suggest  a 
causal  relationship  between  the  delamination  of  mantle  lithosphere  and  the  formation  of 
the  north-south  trending  rift  in  southeastern  Tibet. 

Citation:  Ren,  Y.,  and  Y.  Shen  (2008),  Finite  frequency  tomography  in  southeastern  Tibet:  Evidence  for  the  causal  relationship 
between  mantle  lithosphere  delamination  and  the  north- south  trending  rifts,  J.  Geophys.  Res.,  113 ,  B 103 16,  doi:10.1029/ 
2008JB005615. 


1.  Introduction 

[2]  Since  the  collision  of  India  and  Eurasia  about  50  Ma 
ago,  at  least  2000  km  of  convergence  has  been  accommo¬ 
dated  by  thickening  the  cmst  and  elevating  the  Himalayan- 
Tibetan  Plateau  [Yin  and  Harrison ,  2000].  The  region  has 
been  strongly  deformed  through  a  combination  of  thmst, 
extension,  and  strike-slip  faulting.  Among  the  various  styles 
of  deformation,  north- south  and  northeast- southwest 
trending  rift  zones,  which  indicate  generally  east-west 
extension,  are  enigmatic  features.  They  are  distributed 
mainly  in  the  southern  and  central  Tibetan  Plateau  [Molnar 
and  Tapponnier ,  1978;  Armijo  et  al. ,  1986,  1989;  Rothery 
and  Drury ,  1984],  where  topography  reaches  high  altitudes. 
Field  observations  indicate  that  the  magnitude  of  this 
extension  is  small  (<1%)  [Armijo  et  al ,  1986]  and  rifting 
in  southern  and  central  Tibet  started  about  13.5-14  Ma  ago 
[Yin  et  al ,  1994;  Coleman  and  Hodges ,  1995;  Harrison  et 
al ,  1995;  Blisniuk  et  al ,  2001]. 

[3]  Whether  the  north- south  trending  rifts  are  shallow 
cmstal  features  or  a  surface  manifestation  of  a  coherent 
deformation  throughout  the  entire  cmst  and  mantle  litho¬ 
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sphere  is  a  key  question  in  the  continuing  debate  about  the 
mechanisms  that  form  the  Tibetan  Plateau  [e.g.,  Holt ,  2000; 
Yin ,  2000].  The  late  Cenozoic  east- west  extension  of  the 
plateau  is  commonly  attributed  to  gravitational  collapse  of 
the  thickened  cmst  [Mercier  et  al. ,  1987;  Dewey ,  1988; 
England  and  Houseman ,  1989;  Liu  and  Yang ,  2003],  the 
consequence  of  a  sharp  increase  in  potential  energy  after  an 
abmpt  rise  of  the  plateau  due  to  convective  removal  of  the 
lower  mantle  lithosphere  [England  and  Houseman ,  1989]. 
Numerical  modeling  indicates  that  the  present  high  surface 
elevation  and  small  regional  topographic  slope  within  the 
plateau  require  the  removal  of  the  underlying  mantle  lith¬ 
osphere  [Jimenez-Munt  and  Platt ,  2006].  However,  the  time 
scale  involved  in  the  convective  removal  of  the  mantle 
lithosphere  is  uncertain  because  the  viscosity  is  not  well 
known.  Numerical  simulations  by  different  authors  show 
that  the  duration  of  the  removal  can  vary  drastically  from 
10  Ma  to  a  few  hundreds  of  million  years,  depending  on  the 
rheological  parameters  in  the  simulation  experiments 
[Davaille  and  Jaupart ,  1993;  Conrad ,  2000;  Houseman  et 
al ,  2000;  Morency  et  al ,  2002].  The  roughly  concurrent 
onsets  of  tectonic  processes  in  the  surrounding  regions  have 
led  to  the  suggestion  of  a  rapid  removal  of  the  mantle 
lithosphere  beneath  Tibet  at  about  8  Ma  ago  [Harrison  et 
al ,  1992;  Molnar  et  al ,  1993]. 

[4]  Observations  that  support  the  notion  of  a  vertically 
coherent  deformation  of  the  lithosphere  [Tapponnier  et  al , 
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(a)  Stations  (b)  Events 


Figure  1.  (a)  Topographic  map  of  southeastern  Tibet  with  active  faults  in  the  region  of  the  Himalayan 
Eastern  syntaxis  and  the  locations  of  the  stations  used  in  this  study  (triangles),  (b)  Distribution  of  the 
earthquakes  (circles)  used  in  this  study.  Triangles  represent  the  stations  used.  Numbers  are  the  distance  in 
degrees  to  the  center  of  the  stations. 


1982;  England  and  Houseman ,  1986]  include  the  studies  of 
earthquake  focal  mechanisms  which  show  brittle  east-west 
extension  in  the  mantle  directly  beneath  the  north -south 
trending  rifts  [ Chen  and  Kao ,  1996;  Chen  and  Yang ,  2004] 
and  the  instability  analysis  of  the  spacing  between  the  rifts 
[Yin,  2000].  On  the  other  hand,  the  presence  of  a  low- 
viscosity  mid  and  lower  crust  may  decouple  deformation  in 
the  shallow  crust  from  that  in  the  mantle  lithosphere  [Zhao 
and  Morgan ,  1987;  Bird ,  1991;  Royden ,  1996].  Other 
explanations  for  the  east-west  extension  in  the  Tibetan 
Plateau  include  the  combination  of  the  eastward  extrusion 
of  Tibet  and  strike-slip  faulting  along  the  Himalayan  arc 
caused  by  oblique  convergence  between  India  and  Eurasia 
[Tapponnier  et  al. ,  1981;  Armijo  et  al. ,  1986;  McCaffery 
and  Nabelek ,  1998;  Tapponnier  et  al .,  2001],  and  the 
collisional  stresses  localized  along  the  southern  part  of  the 
Himalayan  arc  [Kapp  and  Guynn ,  2004]. 

[5]  Direct  observational  constraints  on  the  crust  and 
mantle  structure  beneath  the  rifts  are  needed  to  understand 
the  east-west  extension  of  the  Tibetan  Plateau  and  the 
dynamic  processes  that  generate  the  rifts.  In  this  paper,  we 
present  tomographic  images  of  the  lower  crust  and  upper 
mantle  beneath  a  north -south  trending  rift  in  the  region  of 
the  Eastern  Himalayan  Syntaxis,  which  marks  the  eastern 
end  of  the  Himalayan  orogenic  belt  (Figure  la).  The  results 
support  the  notion  of  the  mantle  control  of  the  north- south 
trending  rifts  and  suggest  that  delamination  of  the  mantle 
lithosphere  plays  an  important  role  in  the  rise  of  the  plateau 
in  the  study  area. 

2.  Data 

[6]  The  data  used  in  this  study  come  from  three  sources  in 
southeastern  Tibet:  (1)  the  Namche  Barwa  seismic  experi¬ 
ment,  which  deployed  a  50-station  array  of  broadband 
seismometers  and  a  20-element  array  of  short-period  seis¬ 
mometers  during  2003-2004  [Sol  et  al ,  2007],  (2)  the 
Global  Seismic  Network  (GSN)  station  LSA,  and  (3)  the 


broadband  stations  in  the  Bhutan  experiment,  which 
recorded  concurrently  with  the  Namche  Barwa  array  in 
2003  (Figure  la).  We  used  695  events  occurred  during  the 
operation  period  of  the  stations  in  2003  and  2004,  and  with 
magnitude  greater  than  5.5  (Figure  lb).  With  the  multichan¬ 
nel  waveform  cross  correlation  method  [VanDecar  and 
Crosson ,  1990],  we  measured  differential  travel  times  of 
teleseismic  P  and  S  waves  between  different  stations  in 
short,  intermediate,  and  long  periods:  0.5 -2.0,  0.1 -0.5,  and 
0.03-0.1  Hz  for  P  waves  and,  0.1-0. 5,  0.05-0.1,  and 
0.02-0.05  Hz  for  S  waves.  The  P  differential  traveltimes 
are  measured  on  vertical  component  seismograms  and  those 
of  S  on  transverse  components.  In  order  to  avoid  the 
ambiguity  of  long-period  phase  arrivals  in  the  triplication 
range  and  the  effects  of  the  core-mantle  boundary,  we  have 
considered  only  data  with  epicentral  distance  between  34° 
and  81°. 

[7]  Our  final  data  set  consists  of  ~363 13  P  (19,131, 
10,090  and  7092  for  short,  intermediate,  and  long  periods, 
respectively)  and  ^15043  S  differential  traveltimes  (6629, 
6040,  and  2374  for  short,  intermediate,  and  long  periods, 
respectively),  which  are  then  utilized  to  invert  for  spatial 
variations  in  P  and  S  wave  velocity  perturbations  according 
to  the  3-D  finite  frequency  kernel  formulation  [Dahlen  et 
al ,  2000;  Hung  et  al ,  2004].  For  the  inversion,  we  used 
either  the  combination  of  data  at  all  frequency  bands,  the 
combination  of  only  high-  and  intermediate-frequency  data 
and  only  the  high-frequency  data  to  test  the  reliability  of 
different  data.  The  different  combinations  yield  similar 
results,  though  adding  lower  frequency  data  reduces  the 
variance  reduction  of  the  final  model.  This  is  likely  caused 
by  the  relatively  larger  errors  in  the  long-period  measure¬ 
ments.  The  tomographic  models  discussed  in  this  paper 
come  from  the  inversion  of  the  high-  and  intermediate- 
frequency  data.  Figure  2a  shows  the  average  teleseismic 
differential  traveltime  residuals  at  each  station  for  P  and  S 
waves.  They  are  consistent  with  P  and  S  wave  velocity 
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a)  Mean  traveltime  residuals 
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b)  Station  terms  from  inversion 
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Figure  2.  (a)  Mean  traveltime  residuals  at  each  station  for  (left)  P  wave  and  (right)  S  wave.  Blue 
indicates  an  earlier  than  expected  arrival;  red  indicates  a  later  than  expected  arrival.  Delay  magnitude 
correlates  linearly  with  the  size  of  symbols,  (b)  Station  terms  for  (left)  P  wave  and  (right)  S  wave,  which 
are  solved  simultaneously  in  the  inversion  in  order  to  absorb  travel  time  fluctuations  due  to  uncorrected 
shallow  structures. 


anomalies  obtained  from  the  inversion,  discussed  in 
sections  3-5. 

3.  Methodology 

[8]  Most  of  global  and  regional  tomographic  studies  are 
based  on  ray  theory,  where  seismic  waves  are  assumed  to 
have  an  infinite  frequency  and  the  arrival  time  of  a  body 
wave  phase  depends  only  upon  the  wave  speed  along  the 
geometrical  ray  path  between  the  source  and  receiver. 
However,  because  of  wave  front  healing,  scattering,  and 
other  diffraction  effects,  the  traveltime  of  a  finite  frequency 
seismic  wave  is  sensitive  to  a  three-dimensional  structure 
off  the  ray  path  [Dahlen  et  al .,  2000].  Using  body  wave  ray 
theory  in  conjunction  with  the  Bom  single-scattering  ap¬ 
proximation,  Dahlen  et  al  [2000]  derived  the  formulation 
of  the  Bom-Frechet  kernels  for  a  seismic  phase,  which 
expresses  the  influence  of  velocity  perturbations  upon  a 
finite  frequency  travel  time  shift: 

-///.  ^(x)<5c(x)/c(x)<i3x  (1) 

where  K  is  the  3-D  Frechet  sensitivity  kernel  for  a  shift  6  t , 
measured  by  cross  correlation  of  an  observed  pulse  with  its 


spherical  Earth  synthetics.  It  is  expressed  by  the  formulation 
[Dahlen  et  al ,  2000] 

V-  1  (  R  \ /o°°^3|^»H|2sm(^Ar)^ 

2tt c\crR'R")  JMv.MI  2du,  U 

where  AT  =  T  +  l"  —  T  represent  the  difference  in  travel 
time  for  the  path  with  a  detour  through  a  single  point 
scatterer  between  the  source  and  receiver  (time  I*  +  T") 
relative  to  the  corresponding  direct  ray  path  (time  7);  R ,  R' 
and  R"  are  geometrical  spreading  factors  for  the  unpertur- 
bated  ray,  the  forward  source-to-scatterer  ray  and  the 
backward  receiver-to-scatterer  ray,  respectively.  The  power 
spectmm,  \ssyn\2,  is  for  the  synthetic  seismograms  used  for 
cross  correlation.  Cross-correlating  with  synthetics  in 
different  bandwidths  takes  advantage  of  the  frequency 
dependence  of  the  traveltimes. 

[9]  In  this  study,  we  used  the  finite  frequency  seismic 
tomography,  as  described  by  Hung  et  al  [2004]  and  Yang  et 
al  [2006],  to  retrieve  the  cmstal  and  mantle  stmctures 
beneath  southeastern  Tibet.  Our  measured  data  dt  of  P  and 
S  differential  travel  times  can  be  inverted  to  produce 
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Figure  3.  Trade-off  curve  between  the  model  L2  norm  and 
the  data  variance  reduction  for  different  values  of  the  norm 
damping  parameter  for  P  and  S  waves.  The  larger  square 
and  dot  represent  parameters  used  to  obtain  the  P  and  S 
velocity  models  discussed  in  this  paper. 

tomographic  models  of  velocity  perturbations  using  the 
following  formulation: 

dt  =  [  gi{x)m(x)d3x  (3) 

JD 

where  dh  i  =  1 . .  .N,  represents  the  zth  differential  travel  time 
data;  x  is  the  position  vector  in  3-D  vector  model  space,  D; 
and  gi  represents  the  3-D  Frechet  sensitivity  kernel  relating 
to  the  data  dL  with  the  model  function  m(\). 

[10]  The  crustal  and  mantle  volume  beneath  the  region  of 
study  is  parameterized  with  regular  3D  grids  of  33  x  33  x 
33  centered  at  (93.0°E,  30.0°N),  and  with  a  dimension  of 
18°  in  longitude,  14°  in  latitude,  and  1200  km  in  depth.  This 
results  in  a  grid  spacing  of  ^53  km  in  longitude,  46  km  in 
latitude  and  ~40  km  vertically.  With  this  parameterization, 
equation  (3)  can  thus  be  written  as 

di  =  Gnmi  (4) 

where  dt  is  the  zth  differential  traveltime  data,  Ga  represents 
the  differential  value  of  the  integrated  volumetric  kernels 
contributing  to  the  /th  node,  and  mi  is  the  model  parameter 
at  the  /th  node.  The  inversion  problem  is  resolved  by  the 
standard  damped  least  square  method  [Paige  and  Saunders , 
1982a,  1982b]: 

m  =  (GrG  +  02l) - 1  Grd  (5) 

where  I  is  the  identity  matrix.  The  damping  parameter  6  is 
determined  empirically  through  a  space  of  variance 
reduction  and  the  model  norm  represented  by  a  trade-off 
curve  (Figure  3).  We  choose  the  damping  parameter  that 


yields  an  optimum  variance  reduction  and  a  relatively  small 
model  norm.  The  model  discussed  in  this  paper  is  obtained 
using  a  damping  parameter  that  yields  a  variance  reduction 
of  ~85%  for  P  wave  and  ^80%  for  S  wave. 

[n]  The  traveltime  anomalies  can  be  caused  by  crustal 
thickness,  station  elevations,  and  lateral  velocity  heteroge¬ 
neities.  Because  the  ray  paths  of  teleseismic  body  waves  are 
nearly  vertical  and  crossing  paths  are  rare  at  shallow  depths, 
the  shallow  velocity  structures  are  poorly  constrained  by 
relative  traveltime  data  in  the  inversion.  Therefore,  a  time 
correction  is  needed  in  order  to  reduce  the  tradeoff  between 
crustal  and  mantle  velocity  heterogeneities  in  seismic  to¬ 
mography.  To  correct  for  the  topographic  effects,  we  have 
calculated  the  P  and  S  phase  time  delays  due  to  station 
elevations,  using  P  and  S  ray  parameters  in  the  model 
iasp91  and  a  velocity  of  5.8  km/s  and  3.36  km/s  for  P 
and  S  waves,  respectively.  The  measured  differential  trav- 
eltimes  are  then  subtracted  by  the  calculated  P  and  S  time 
delays  for  the  topographic  correction.  To  account  for  the 
thick  crust  beneath  Tibet,  we  first  replace  the  crust  in  the 
iasp91  model  with  the  average  crust  structure  beneath 
the  stations  from  the  crustal  model  CRUST2.0  [Bassin  et 
al. ,  2000].  Then,  an  additional  free  term  at  each  station  is 
incorporated  into  the  inversion  to  absorb  traveltime  shifts 
caused  by  other  shallow  heterogeneity.  Figure  2b  shows  that 
the  simultaneously  solved  station  terms  yield  a  consistent 
pattern  among  P  and  S  wave  velocity  models.  The  maxi¬ 
mum  difference  between  the  station  terms  is  ~0.36  s  and 
^0.42  s  for  P  and  S  models,  respectively. 

4.  Resolution  Tests 

[12]  We  performed  different  types  of  resolution  tests  to 
evaluate  our  data  coverage  and  the  ability  of  our  inversion 
technique  to  recover  the  mantle  structure.  For  this  purpose, 
the  synthetic  traveltimes  are  computed  by  multiplying  the  G 
matrix  (equation  (4))  with  different  input  velocity  models: 
A tsyn  =  G  •  A csynth.  The  inversion  was  then  performed  using 
the  same  damping  parameter  as  in  the  inversion  of  real  data. 
The  resolution  is  considered  good  when  input  structures  are 
well  recovered,  as  illustrated  in  Figures  4-6.  In  the  first  set 
of  the  tests  concentrating  on  the  horizontal  resolution,  we 
conducted  the  checkerboard  resolution  tests  using  succes¬ 
sively  smaller-scale  input  anomalies  with  horizontal  sizes  of 
200  km  x  200  km,  150  km  x  150  km,  and  100  km  x 
100  km.  In  each  case,  the  input  models  are  composed  by  a 
velocity  perturbation  of  ±2%  for  P  wave  and  ±4%  for  S 
wave,  which  alternate  horizontally  (Figure  4a  for  P  wave 
and  Figure  5a  for  S  wave).  We  estimated  from  the  different 
tests  that  the  resolution  is  good  down  to  ^450  km  depth  for 
anomalies  greater  than  and  equal  to  150  km,  and  ^300  km 
for  100-km-sized  anomalies.  The  horizontal  dimensions  of 
the  interpreted  features  in  our  tomographic  model  are  about 
150  km  or  more.  In  the  second  set  of  the  tests  for  both 
horizontal  and  vertical  resolutions,  the  input  anomalies 
alternate  both  horizontally  and  vertically.  These  tests  show 
that  the  minimum  size  of  the  anomaly  which  can  be 
recovered  vertically  is  ~120  km  in  the  upper  400-km 
mantle  (Figure  4b  for  P  wave  and  Figure  5b  for  S  wave). 

[13]  In  order  to  obtain  a  quantitative  assessment  of  the 
resolution  tests,  we  calculated  the  correlation  coefficient 
between  the  input  and  output  models  within  the  well 
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Figure  4.  (a)  Checkerboard  resolution  tests:  (a)  Horizontal  slices  of  the  input  and  retrieved  P  wave 
velocity  models  at  depths  112  km  and  338  km.  The  magnitudes  of  the  input  anomalies  are  ±2.0%;  the 
sizes  of  the  anomalies  are  (left)  MOO  km  x  200  km,  (middle)  ~150  km  x  150  km  and  (right) 
~100  km  x  100  km.  The  rectangle  with  dashed  lines  on  the  horizontal  slice  at  112  km  depth  marks  the 
area  which  we  consider  well  resolved  and  used  for  Figure  6.  (b)  (left)  Vertical  slice  along  the  path  A-B 
through  the  input  and  retrieved  P  wave  velocity  model.  The  horizontal  size  of  the  anomalies  is 
MOO  km  x  200  km  and  the  vertical  size  is  M20  km.  (right)  The  same  vertical  slice  for  the  input 
anomalies  with  the  horizontal  scale  of  M50  km  x  150  km  and  the  vertical  scale  of  M20  km. 


resolved  area  (rectangle  with  dashed  lines  in  Figure  4a)  as  a 
function  of  depth  (Figure  6a),  as  well  as  the  standard 
deviation  of  the  difference  between  the  input  and  output 
models  as  a  function  of  depth  (Figure  6b),  for  the  three 


cases  of  the  input  velocity  anomalies  above.  For  the  P  and  S 
models  with  the  two  larger-sized  anomalies,  the  correlation 
coefficients  between  the  input  and  ouput  models  are  very 
good  in  the  upper  450  km  depth  (>0.9).  The  standard 
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Figure  5.  Checkerboard  resolution  tests:  (a)  Horizontal  slices  of  the  input  and  retrieved  S  wave  velocity 
models  at  depths  112  km  and  338  km.  The  magnitudes  of  the  input  anomalies  are  ±.0%;  the  sizes  of  the 
anomalies  are  (left)  ^200  km  x  200  km,  (middle)  ~  150  km  x  150  km  and  (right)  ^100  km  x  100  km. 
The  rectangle  with  dashed  lines  on  the  horizontal  slice  at  1 12  km  depth  marks  the  area  which  we  consider 
well  resolved  and  used  for  Figure  6.  (b)  (left)  Vertical  slice  along  the  path  A-B  through  the  input  and 
retrieved  S  wave  velocity  model.  The  horizontal  size  of  the  anomalies  is  ^200  km  x  200  km  and  the 
vertical  size  is  ~120  km.  (right)  The  same  vertical  slice  for  the  input  anomalies  with  the  horizontal  scale 
of  ~  150  km  x  150  km  and  the  vertical  scale  of  ~120  km. 
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Figure  6.  (a)  Correlation  coefficient  between  the  input  and  output  (left)  P  and  (right)  S  models  as  a 
function  of  depth  for  the  cases  in  which  the  sizes  of  the  input  anomalies  are  ^200  km  x  200  km  (line 
with  squares),  ~  1 50  km  x  150  km  (line  with  stars),  and  ~100  km  x  100  km  (line  with  triangles).  The 
calculations  are  for  the  area  that  we  consider  well  resolved  (rectangle  with  dashed  lines  in  Figure  4a). 
(b)  Standard  deviation  of  the  difference  between  the  input  and  output  models  as  a  function  of  depth  for 
the  three  cases  of  the  input  anomalies. 
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Figure  7.  (a)  P  and  S  wave  velocity  perturbations  ( dlnVP  and  dlnVs)  relative  to  the  reference  model 
(CRUST2.0  and  iasp91)  at  different  depths:  75  km,  112  km,  188  km,  and  375  km.  The  rectangle  with 
dashed  lines  marks  the  area  which  we  consider  well  resolved  (also  see  Figures  4a  and  5a).  (b)  Cross 
section  A-B  through  P  and  S  wave  tomographic  models,  along  the  Indus-Yarlung  suture  and  across  the 
rift  in  southeastern  Tibet. 


deviations  of  the  difference  between  the  input  and  output 
models  are  respectively  0.7%  and  1.3%  for  the  P  and  S 
models.  This  means  that  about  65%  of  the  amplitudes  of  the 
input  models  at  these  scales  can  be  recovered.  For  the  case 
in  which  the  size  of  the  input  anomalies  is  ~100  km  x 
100  km,  the  recovery  in  amplitude  is  about  55%  in  the 
upper  300  km  and  degrades  gradually  at  greater  depths. 

5.  Results  and  Discussions 

[14]  Figure  7a  shows  P  and  S  wave  velocity  structures  in 
the  crust  and  upper  mantle  beneath  southeastern  Tibet  at 
four  different  depths.  There  is  a  prominent  north- south 
trending,  low-velocity  structure  at  ~92°  longitude  from  the 
lower  crust  to  ^300  km  depth.  This  structure  is  observed  on 
both  the  P  and  S  velocity  models.  It  extends  across  the 
Indus-Yarlung  Suture  and  coincides  strikingly  well  with  a 
rift  on  the  surface.  The  shallow  (75  km  and  112  km) 
velocity  structure  (Figure  7a)  is  in  agreement  with  Pn  wave 


tomography  [Liang  and  Song ,  2006],  which  also  shows  a 
north -south  trending,  low- velocity  anomaly  in  the  same 
region,  albeit  at  a  coarse  resolution  and  centered  somewhat 
west  of  the  low-velocity  feature  in  this  study.  At  depth 
above  250  km,  there  is  a  high-velocity  anomaly  east  of  the 
north -south  trending,  low- velocity  feature  and  north  of  the 
surface  trace  of  the  main  boundary  thrust  fault  (MBTF) 
between  the  Indian  and  Eurasian  plates.  We  associate  this 
high-velocity  anomaly  to  the  Indian  lithosphere  underthrust¬ 
ing  beneath  the  plateau.  While  we  agree  with  the  interpre¬ 
tation  of  Pn  tomography  that  the  Indian  lithosphere 
advances  further  north  near  the  eastern  comer  than  to  the 
west  [Liang  and  Song ,  2006],  the  high-velocity  anomaly 
associated  with  the  Indian  lithosphere  appears  to  be 
bounded  by  the  Jiali  fault.  There  is  no  evidence  for  a 
subducted  Indian  lithosphere  at  depths  greater  than  250  km 
beneath  the  study  area. 

[15]  West  of  the  north- south  trending,  low- velocity 
anomaly  is  a  high-dipping-angle,  high-velocity  anomaly 
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Figure  8.  Resolution  tests  using  tabular  structures.  The  cross  section  A-B  is  the  same  as  that  in  the 
models  presented  in  Figure  7b.  (top)  Input  P  and  S  wave  velocity  models  (+2%  and  +4%,  respectively) 
and  the  corresponding  Vp/Vs  ratio  (—2%).  The  VP/VS  ratio  is  obtained  by  adding  the  VP  and  Vs 
perturbations  from  the  inversions  to  the  P  and  S  reference  models,  dividing  the  perturbed  P  and  S 
velocities,  subtracting  the  reference  VP/Vs,  and  then  calculating  the  VP/VS  change  in  percentage,  (bottom) 
Retrieved  P  and  S  wave  velocity  models  and  the  corresponding  VP/VS  ratio. 


that  extends  into  the  upper  mantle  transition  zone,  as  shown 
on  the  cross  sections  in  Figure  7b.  Resolution  tests  using 
synthetic  slabs  with  a  thickness  of  ^120  km,  placed  either 
in  the  upper  and  lower  parts  of  the  upper  mantle,  show  little 
vertical  smearing  in  the  recovered  structures  (Figure  8).  We 
conclude  that  the  observed  high-velocity  anomaly  is  a 
robust  feature.  This  anomalous  structure  is  also  more  or 
less  tabular  as  opposed  to  cylindrical.  The  possible  causes 
for  seismic  velocity  anomalies  in  the  mantle  are  variations 
in  temperature,  chemical  composition,  and  volatile  contents. 
Thermal  variations  alone  cannot  fully  account  for  the 
observed  high-velocity  anomaly,  which  coincides  with  a 
low  Vp/Vs  ratio  (Figure  7b).  A  temperature  reduction  of 
400°  C  for  a  wide  range  of  upper  mantle  compositions 
causes  approximately  a  —0.5%  change  in  VP/VS  [Hacker 
and  Abers,  2004;  Boyd  et  al .,  2004],  much  smaller  than  the 
observed  anomaly  (—1.5%).  Resolution  tests  show  that  for 
features  comparable  to  the  interpreted,  delaminated  mantle 
lithosphere  and  the  Indian  mantle  lithosphere  beneath  the 
eastern  indentor,  the  inversion  yields  well-recovered  Vp/Vs 
ratios,  with  magnitudes  slightly  less  than  those  of  the  input 
anomalies  (Figure  8).  So  the  magnitude  of  the  real  Vp/Vs 
ratio  of  the  high-angle,  high-velocity  anomaly  is  likely 
greater  than  the  observed  value  in  Figure  7b,  requiring  a 
compositional  change,  which  could  be  accounted  for  by  a 
refractory  mantle  depleted  of  volatiles  [Hacker  and  Abers , 
2004;  Boyd  et  al ,  2004].  The  fact  that  the  high-velocity 
anomaly  in  the  shallow  mantle  beneath  the  eastern  indentor, 
which  we  identify  as  the  Indian  mantle  lithosphere,  also  has 
a  low  Vp/Vs  of  a  similar  magnitude  supports  the  interpre¬ 
tation  that  the  high-dipping-angle,  high-velocity  anomaly  is 


a  sunken  mantle  lithosphere.  The  north -south  orientation  of 
this  feature  excludes  the  possibility  that  it  is  a  detached 
Indian  plate.  We  suggest  that  the  low-velocity  anomaly 
immediately  east  of  and  above  this  high-velocity  feature 
represents  the  mantle  asthenosphere  that  filled  the  void  left 
by  the  sunken  Eurasian  mantle  lithosphere  (Figure  9).  The 
tabular  nature  of  the  high-angle,  high-velocity  anomaly,  the 
asymmetry  of  the  low-velocity  anomaly  above  it  and  the  rise 
of  this  low-velocity  anomaly  to  the  base  of  the  crust  at 
~70  km  suggest  a  process  in  which  the  entire  mantle 
lithosphere  peeled  off  as  opposed  to  a  partial  removal  of 
the  mantle  from  the  base  of  a  thickened  lithosphere  by 
viscous  flow.  The  shallow  asthenosphere  and  associated 
heating  and  weakening  of  the  overlying  crust  could  thus 
initiate  and  localize  the  rifts  observed  on  the  surface  of  the 
Tibetan  Plateau,  at  least  in  the  region  of  the  Eastern 
Himalayan  Syntaxis. 

[i6]  Roughly  coincident  with  the  onset  of  the  east- west 
extension,  potassic  volcanism  became  widespread  on  the 
Tibetan  Plateau.  In  southern  Tibet,  volcanism  may  have 
started  slightly  earlier  (^16-20  Ma)  and  then  ceased 
around  10  Ma  [WUliams  et  al .,  2004].  The  small  magnitude 
of  the  east-west  extension  is  insufficient  to  cause  decom- 
pressional  melting  of  the  asthenosphere  [McKenzie  and 
Bickle ,  1988].  Geochemical  modeling  shows  that  the  mag¬ 
mas  in  southern  Tibet  were  derived  from  a  small  degree 
(<2%)  of  partial  melting  of  metasomatized  (phlogopite) 
peridotite  in  the  spinel  stability  field  (at  depths  of  65- 
80  km)  [Williams  et  al ,  2004].  Our  tomographic  results  are 
consistent  with  the  geochemical  analysis.  The  delamination 
of  the  mantle  lithosphere  brought  the  subcontinental  mantle 
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Figure  9.  Three-dimensional  image  of  the  P  wave  velocity  model  with  our  interpretation  of  the 
observed  structures.  The  topography  is  not  to  scale. 


lithosphere  into  direct  contact  with  astheno spheric  temper¬ 
atures,  so  that  metasomatized  peridotite,  which  has  a  lower 
solidus,  underwent  partial  melting  right  beneath  the  crust  or 
released  its  volatiles  into  the  shallower  mantle.  The  tabular 
geometry  of  the  north- south  trending  low- velocity  anomaly 
suggests  that  the  release  of  volatiles  from  the  sunken  mantle 
lithosphere  may  contribute  significantly  to  the  velocity 
reduction  above  it  since  a  buoyant  and  low-viscosity  up- 
welling  tends  to  form  a  cylindrical  not  tabular  anomaly. 
Melting  would  cease  or  be  greatly  reduced  once  the  vola¬ 
tiles  in  the  peridotite  are  depleted  by  the  melting  process. 
The  fact  that  extension  postdates  the  earliest  magmatism  by 
a  few  million  years  [  Williams  et  al .,  2004]  further  suggests  a 
causal  relationship  between  the  mantle  processes  and  the 
surface  rifts. 

[17]  The  thickening  and  deformation  of  the  lithosphere 
caused  by  the  Indo-Eurasia  collision  is  certainly  a  likely 
geological  setting  for  the  initiation  of  the  lithospheric 
delamination  [Bird,  1979],  though  the  exact  mechanism 
for  the  initiation  of  the  delamination  of  the  mantle  litho¬ 
sphere  is  unclear.  Shear  heating  through  viscous  dissipation 
in  the  lithosphere  likely  plays  a  role  in  weakening  the 
lithosphere  [Kincaid  and  Silver ,  1996;  Schott  et  al ,  2000] 
and  the  initiation  of  the  delamination  at  a  localized  shear 
zone.  Numerical  simulations  on  the  dynamics  of  the  mantle 
lithosphere  show  that  mantle  delamination  and  detachment 
likely  occur  if  the  lithosphere  is  substantially  thickened  and 
that  the  process  is  strongly  controlled  by  the  rheology  of  the 
lower  crust  [Schott  and  Schmeling,  1998;  Morency  and 
Doin ,  2004].  These  experiments  show  that  a  low- viscosity 
and  hot  lower  crust  favors  especially  the  delamination 
phenomenon  in  a  thickened  lithosphere  and  the  Tibetan 
Plateau  fulfils  well  these  conditions.  Several  recent  Hima- 
layan-Tibetan  tectonic  investigations  have  proposed  a  weak 


middle-lower  crust,  which  allows  crustal  flow,  to  explain  the 
building  of  the  southeastern  margin  of  the  Tibetan  Plateau 
[Clark  and  Royden,  2000;  Beaumont  et  al .,  2001; 
Schoenbohm  et  al ,  2006].  This  idea  is  supported  by  results 
from  geophysical  studies  that  show  evidence  for  localized 
partial  melt  within  the  middle  crust  in  Tibet  [Kind  et  al , 
1996;  Wei  et  al ,  2001]  and  a  radial  anisotropy  that  can  be 
caused  by  channel  flow  within  the  mid-to-lower  Tibetan 
crust  [Shapiro  et  al ,  2004].  However,  contrary  to  the  notion 
that  the  north -south  rift  zones  in  the  Tibetan  Plateau  are 
shallow  features,  formed  by  the  eastward  motion  of  the 
shallow  crust  that  are  decoupled  from  the  mantle  lithosphere 
by  a  low-viscosity  lower  crust,  our  results  suggest  an  upper 
mantle  origin  of  the  rift  zones,  at  least  in  southeastern  Tibet, 
where  mantle  lithosphere  delamination  plays  a  key  role  in 
the  process  of  the  rise  of  the  plateau. 

6.  Conclusion 

[is]  Finite  frequency  tomography  of  teleseismic  P  and  S 
traveltimes  shows  the  presence  of  a  low- velocity  anomaly  in 
the  crust  and  upper  mantle  down  to  ~300  km  depth  beneath 
a  north- south  trending  rift  zone  in  southeastern  Tibet.  This 
low-velocity  anomaly  is  situated  above  a  tabular,  high- 
dipping-angle,  high-velocity  anomaly  that  extends  into  the 
upper  mantle  transition  zone.  The  VP/VS  ratio  of  this  high- 
velocity  anomaly  suggests  that  temperature  variations  are 
not  the  only  cause  and  a  highly  melt-depleted  mantle  is 
required.  These  results  provide  evidence  for  the  mantle 
lithosphere  delamination  and  its  link  to  north -south  trend¬ 
ing  rifts  in  southeastern  Tibet.  Studies  on  an  extended 
region  in  southern  Tibet  are  thus  needed  to  understand 
whether  this  phenomenon  is  limited  to  the  study  area  or 
occurs  in  a  broad  region  where  mantle  lithosphere  delam- 
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ination  plays  an  important  role  in  the  elevation  and  defor¬ 
mation  of  the  Tibetan  Plateau. 
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